{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tails in social data\n",
    "Fiona Pigott\n",
    "\n",
    "This is going to be a multi-part tutorial on dealing with long-tail distributions in a social data context. This first tutorial will define some vocabulary, look a few different distributions, glance at some empirical data, and introduce a nifty little python package [powerlaw](https://github.com/jeffalstott/powerlaw/) for fitting distributions (power law, lognormal and exponential) by Jeff Alstott. \n",
    "\n",
    "## Resources\n",
    "- [Not every long tail is power law!](http://www.philippsinger.info/?p=247) by Philipp Singer of the [GESIS - Leibniz Institute for the Social Sciences](http://www.gesis.org/home/) \n",
    "- [Power-Law Distributions in Empirical Data](http://dx.doi.org/10.1137/070710111) by Aaron Clauset, Cosma Rohilla Shalizi, and M. E. J. Newman\n",
    "- [Slides](http://users.cms.caltech.edu/~adamw/heavytails.html) and [blog posts](https://rigorandrelevance.wordpress.com/2013/11/15/a-tale-of-heavy-tails/) on _The Fundamentals of Heavy Tails: Properties, Emergence, and Identification_ by  Adam Wierman, Bert Zwart, JK Nair of Caltech\n",
    "- _An Introduction to Heavy-Tailed and Subexponential Distributions_ by Sergey Foss, Dmitry Korshunov, Stan Zachary\n",
    "- And of course, [Wikipedia](https://www.wikipedia.org/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plotting library\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# mathematics\n",
    "from math import exp, pi, sqrt, log\n",
    "from numpy import linspace, hstack\n",
    "from collections import Counter\n",
    "from random import sample\n",
    "from random import uniform\n",
    "from functools import reduce\n",
    "\n",
    "# special powerlaw fitting library\n",
    "import powerlaw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# x values\n",
    "x = linspace(0,5,1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heavy tails\n",
    "\n",
    "Confession: I absolutely did not know this until I read the Wikipedia article\n",
    "\n",
    "**Definition: A non-negative function (usually tending to zero) is heavy-tailed\n",
    "if it fails to be bounded by a decreasing exponential function.**\n",
    "\n",
    "There are many subtypes of heavy tailed distributions, **one** of which is a long tailed distribution. If you're like me, you just now learned that \"long\" and \"heavy\" are not interchangeable when referring to tails.  ¯\\\\\\_(ツ)_/¯"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential Distribution\n",
    "The exponential distribution is often used to model the time between events (say, the time between Tweets from a particualr user). This is only a valid model if:\n",
    " - Events (Tweets) occur independently.\n",
    " - The rate at which events occur is constant. \n",
    "\n",
    "The PDF of the exponential distribution is: $ f(x) = \\lambda e^{-\\lambda x} $, where $\\lambda$ is commonly referred to as the rate parameter, or the average amount of time between two events (Tweets).\n",
    "\n",
    "The exponential distribution is **not** a long-tailed distribution, but any function with more mass in the tail (tail = $x \\rightarrow \\infty$) than the exponential is a **_heavy tailed_** distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAGHCAYAAAAdnkAlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xm8XHV9//HXB0QWkaggi4KYoEX8aZFEEVxuFBSqom0V\nDWktVhSKxUBDbUVajYBLKwVaUAStstQaoNYFcUEBJVqFSCJQFdwSBDVsotGyKMvn98c5Y04mMzcz\n587cmTvzej4e88jMmXPOfO5J4L7nu53ITCRJknptk0EXIEmSRpMhQ5Ik9YUhQ5Ik9YUhQ5Ik9YUh\nQ5Ik9YUhQ5Ik9YUhQ5Ik9YUhQ5Ik9YUhQ5Ik9YUhQ1JtEfHViPhKzWMfioh31Dx2fnn8RGXbuRGx\nus75anz+TRHx0crr15X1zJ2mz6993aXpZMjQ2Kv8gmj1eDAi9h50jYMUEXtExJKIeEKLtxN4aLpr\nqnx28+uuatnIzzaZh9p8fs8M8XWXOvawQRcgDYkE3g7c1OK9H01vKUPnqcAS4CvAzU3vvXj6y2nr\njXT/xWmyn20yu9P/X/Iz5bpLbRkypHW+mJkrB13EEArafEvPzAemuZa2MvNB4MEuD2v7s7XcOWKL\nzLwvM+/v8nPqmBHXXZqM3SVShyLihLL75IVN2z8cEb+NiKeXrxvjBV4TEe+JiDUR8X8R8ZmI2LnF\neV8dEddExD0RcUdE/EdEPK5pn3Mj4jcR8biI+HT5/PaIODkiomnfiIi/iYjvRMS9EXFrRJwVEY9q\n2u+miLg4Ip4bEVeX+/44Iv6iss/rgIvKl1+tdCFNlO9/NSKuqOy/WUScWP48vyp/7mUR8YJaF704\n5+PLn/n/IuK2iDgV2Jzil3DzNVrdtO2QspZfR8TaiLg+IhZ1+LM1rs8BEfGtiLgPOKLy3kfZ0CMi\n4uyIuLP8vPNaXPeWY1Gq5+z2upfbHhsRHyn/vu+NiGsj4tCmfXYtz3VsRBweET+KiPsiYnlEPLPd\n34FUly0Z0jqzImLbpm2ZmXeVz08CDgI+EhFPz8y7I+JA4A3AP2Tm/zYd+w8UTer/BGwPLAa+HBHP\nyMzfAkTEXwIfBa4GjgN2AP4GeE5E7JWZv27UQfGl4FLgKuBvgRcBx1J055xd+dwPAYeW5/03YDaw\nCHhGRDy3/MbfOOeTgf8CPgKcCxwGnBMR12TmDcAy4PTy+HcBN5bH3lA5R9U25TmWlnU8srw+X4yI\nvTPzeroQEVsAVwA7lz/LGuAvgP1afHZWt0XEi4GPA18G/r7cvAfwHOCMDn+2p5TnOLv8eb7f5ueG\nIvS8H/glRTfHHwBHAU8AXthi/2bVc3Z13cvr9FVgt/Jnuwl4NXBuRMzKzDOaPuvPga2Bs8pzvRX4\n74iYU/n3IU1dZvrwMdYP4HUUYaDV456mff8fcB/FL51ZwE8pfulvUtlnfnnszcBWle0Hl9vfXL5+\nGHArcC3w8Mp+Ly33W1LZdg5FV8DxTfWsAJZXXj+vPHZB034vLrcfUtm2ujzncyrbtgPuBd5X2faq\ncr+JFtfuK8AVldcBPKxpn20owsGHm7Y/BLxjI383x5Sf/crKti2AHzTXVF6jVZXXpwF3beT8k/1s\njevzojbvfbTFv6GrgU0r299SnuOgjf3cLc7ZzXVvXKfq3++mwP8Aa4FHlNt2LT//dmCbyr4vL49/\n6aD/e/QxWg+7S6RCAm+iaB2oPl6y3k6Z36X4lno4RavCY4DXZWarQYDnZeY9lWM/QfHL9qXlpmdR\ntHCcmZm/q+z3eYpvri9rcc6zm15/DZhTeX0w8Cvg8ojYtvEAvg38Hxt+o/5eZn6j8tl3Unxbn0MN\nWXgAft9t82jg4cA1QJ3pnS8B1mTmJyufcR9Fq8LG/ArYumxtqmt1Zl7Wxf4fyvVbAj5I+ct7CjV0\n4iXArZl5QWNDWcfpFC0W85v2vyDXtZJB8e8oqPn3LrVjd4m0zreys4GfJwOHUISE4zPz+232azUr\n5UcU3yahaEZPim/lzW4Entu07b7M/EXTtl8Cj668fjLwKIpvqs2SItRUtZpR0XzOrpTjCY6l6GrY\nrPLWqhqn25XW17HdNa86k6LL4PMR8XPgS8BFmXlpF5/fzbobSVOtWXSprWHd33m/7Ar8sMX2GyjC\nQ/Pn31J9kZm/Kof21P57l1oxZEjd243ilznA07s8Nto870QnfeWbALcBf9bm/Hd0eM5uaysOingt\nRbfFJ4H3UYSdB4Hjqfctud0Mi43Wl5l3RMQzgAMpvum/BHh9RJyXma/v8PPv7bjS9jq9lptOw2c0\n9PTvXWrHkCF1oZzJcS5FP/dpwD9ExCcy89Mtdn9yi227AdeVz2+i+J/67hSD9qp2B35So8QfA/sD\n38hycGkPdLPI1KuAH2fmwdWNEXFizc++CXhai+27d3Jw2XXzufJBRHwQOCIiTsrMVfR2Aa2g+Du/\n8vcbIh4B7AhcUtnvlxStTVT22wzYqbn8Lj77JloH3j3KP+v8W5KmzDEZUnf+FtiHYkzGOygG1n0w\nIh7TYt9DI2LrxouIeDXFL5LPl5uuofimf2T5S6ax30sofjlUfzF16iKKLw+tpkhuGhGzapzzbopf\noI/a2I4U35CbZz48G9i3xudCca12iohXVc63FcX1n1Sbv5PGDKDNyz+7+dk6cUREVL+8/TVFC8Xn\nK9t+DEywviPZsCWjm9o+D+wYEQsaGyJiU4rZKb+hEnyk6WRLhlQI4KURsUeL976RmavL904EzikH\nZxIRr6eYHfJBYEHTcXcBX4+Icyi+zR5DMf7i36H4lh0Rb6WYarosIpaW+x1NMX7hX7v9ITJzWUSc\nDRxXdhV8CbifYjrlweW5PznJKVq5liI8vLVc8+G3wOXlINFmlwCvjIhPU7QezAH+CvguxQDEbn0Y\neDPwH+U6Do0prHd3cOy/l0HjCopZQE8sz3VtFtNzu/3ZOvFwikG3F1GMSXkT8LXMrAbGfwfOiohP\nUEyv3RM4gA27srqp7UMU1/nc8jrdRDEeZV/gmMzs5HpJvTfo6S0+fAz6QTH98MFJHodStPpdTfE/\n70c2Hb+o3O/g8vX88vVrKNY4WEMxs+MzwM4tPv9gilaNeyh+0ZwH7NS0zznA2hbHLgEeaLH9DcDy\n8nN/RfEL6z3ADpV9VgGfaXHsVyh+mVW3HUYxsPB3VKZVttn3reW57yl/rpeU9f+4ab8Hgbd38Pez\nM/Apim/ktwGnUEzJbTWF9ceV138KfKG8/vdSDOL8ALB9hz/b6lbXp3LtPtLi39DzKALnnRRdaucB\nj2o6Nsq/i9vKn+lzFGuZrHfOGtd9O4oAc1v5814L/EXTPruW51nc4mfq6O/Dh49uHpHZyy5JSREx\nn+KXwMFZmXopSeNmKMZkRMTzy+V7f1YuefuKDo55QUSsKJfE/UE5bU6SJA2JoQgZwCMomvaOooMR\n1RHxRIq+38sp+jP/jaL/1TsTSpI0JIZi4GdmfhH4Ivx+iuDGvIli+eDG/Qi+HxHPo7w3RH+qlLpi\nP6SksTcUIaOGfYDmpX4vpVi3QBqozLySqS2sJEkjYVi6S7q1I8UI6qrbgG0iYvMW+0uSpGk2U1sy\nWml0s7Rspi5vEnUgxRTE+6apJkmSRsEWFGvNXJob3kOprZkaMm4Fdmjatj3w66zczbLJgcB/9rUq\nSZJG258DH+9055kaMr5J0y24KVbM++Ykx9wEcNppH2NiotWijuqHxYsXc9ppDpWZTl7z6ec1n35e\n8+l1ww038NrXvhbK36WdGoqQUd5E6Ems6/KYExF7Andl5i0R8V7gcZnZWAvjLODNEfHPFEsy70+x\nauJLJ/mY+wAmJvZg7ty5/fgx1MKsWbO83tPMaz79vObTz2s+MF0NNxiWgZ/PBL4NrKAYU3EKsBI4\noXx/R2CXxs6ZeRPwMuBFFOtrLAbekJnNM042sHJlL8uWJEntDEVLRjnlr23gyczXtzlmXrefdfjh\nsP/+MHt2t0dKkqRuDEtLxrT58IcNGJIkTYexCxmHHw7Llg26ivGxcOHCQZcwdrzm089rPv285jPD\n2NyFNSLmAis+/OEVHH74XFatskVDkqROrFy5knnz5gHMy8yORzeOXUvG3Llw5ZUGDEmS+m3sQsbP\nfgbz59tlIklSv41dyHj844uWjPnzYfXqQVcjSdLoGruQATAxYZeJJEn9NpYhY/Vqu0wkSeq3sQwZ\ns2fbZSJJUr+NZcgAu0wkSeq3sQ0ZdplIktRfYxsy7DKRJKm/xjZkgF0mkiT101iHDLtMJEnqn7EO\nGXaZSJLUP2MdMgB22WXQFUiSNJrGPmQ0WjMclyFJUm+NfchwXIYkSf0x9iHDcRmSJPXH2IcMcCqr\nJEn9YMjALhNJkvrBkIFdJpIk9YMho2SXiSRJvWXIKNllIklSbxkySnaZSJLUW4aMClf/lCSpdwwZ\nFbNnw6pVg65CkqTRYMhoYc4cu0wkSZoqQ4YkSeoLQ0YTu0wkSeoNQ0YbdplIkjQ1hoxJ3HLLoCuQ\nJGnmMmS04JoZkiRNnSGjDZcZlyRpagwZbbjMuCRJU2PIaMMuE0mSpsaQMQmXGZckqT5DxiQarRmO\ny5AkqXuGjEk4LkOSpPoMGZNwXIYkSfUZMjbCcRmSJNVjyNgI72UiSVI9howOeS8TSZK6Y8jogvcy\nkSSpc4aMDjgAVJKk7hkyOuS9TCRJ6o4ho0OumSFJUncMGR2yy0SSpO4YMrrQWDPDAaCSJG2cIaML\ntmZIktQ5Q0aXHAAqSVJnhiZkRMRREbE6Iu6NiKsi4lkb2f9vIuLGiLgnIm6OiFMjYvN+1+kAUEmS\nOjMUISMiFgCnAEuAvYDrgEsjYrs2+/8Z8N5y/6cAhwELgHf3u1a7TCRJ6sxQhAxgMXB2Zp6fmTcC\nRwL3UISHVvYFvp6ZF2bmzZl5GbAU2Hs6ivWmaZIkbdzAQ0ZEbAbMAy5vbMvMBC6jCBOtfAOY1+hS\niYg5wEuBz/W32oI3TZMkaeMeNugCgO2ATYHbmrbfBuze6oDMXFp2pXw9IqI8/qzM/Oe+Vtpkzpwi\nbDgIVJKkDQ28JWMSAWTLNyJeABxP0a2yF/BK4KCI+Mdpq06SJE1qGFoy7gQeBHZo2r49G7ZuNJwI\nnJ+Z55SvvxsRWwNnA++a7MMWL17MrFmz1tu2cOFCFi5c2G3dQLEwly0ZkqRRsXTpUpYuXbretrVr\n19Y6VxTDHwYrIq4Crs7MY8rXAdwMnJ6ZJ7fY/xrgy5n5tsq2hcC/A1tnix8qIuYCK1asWMHcuXN7\nUveyZcUsE7tMJEmjbOXKlcybNw9gXmau7PS4YWjJADgVOC8iVgDLKWabbAWcCxAR5wM/zczjy/0/\nCyyOiGuBq4EnU7RufKZVwOgXZ5lIktTeUISMzLyoHMh5IkW3ybXAgZl5R7nLzsADlUNOAh4q/3w8\ncAdwMTCtYzKcZSJJUntDETIAMvNM4Mw27+3X9LoRME6ahtI2ylkmkiRtaJhnl8wo3plVkqT1GTKm\nyGXGJUlqzZDRAw4AlSRpQ4aMHnAAqCRJGzJk9NCcOXaZSJLUYMjoMQeASpJUMGT0iANAJUlanyGj\nhxoDQG3NkCTJkNFTtmZIkrSOIaPHnM4qSVLBkNFjjdYMlxiXJI07Q0aPrV5ddJcsWzboSiRJGixD\nRo85LkOSpIIhow+cZSJJkiGjL2zNkCTJkNE3tmZIksadIaNPbM2QJI07Q0YfuWaGJGmcPazbAyJi\nc+DZwK7AVsAdwLcz0+/rTbwFvCRpnHUcMiLiucAxwMuBzYC1wL3AY4DNI2IV8CHgrMz8TR9qnbHm\nzCnChgt0SZLGSUfdJRHxGeBC4CbgAOCRmbltZu6cmVsBTwbeBewP/CAiXtynemcsB4BKksZNpy0Z\nXwAOzsz7W72ZmauAVcB5EbEH8Lge1TfjVQeA2pohSRonHbVkZOZZ7QJGi31vyMzLp1bWaHE6qyRp\nHHU9uyQiXjjJe381tXJGk9NZJUnjqM4U1i9GxMkRsVljQ0Q8NiI+C/xT70obLbZmSJLGTZ2Q8QLg\nT4FvRcRTI+JlwHeAbYBn9LC2kWJrhiRp3HQdMjLzm8CeFMFiJfAp4FTgBZn5k96WN1pcnEuSNE7q\nrvi5O/BM4KfAA+XrrXpV1Kizy0SSNA7qDPw8Dvgm8GXgacCzgL2A6yNi396WN1rsMpEkjZM6LRnH\nAH+SmYsy877M/C6wN/BJ4Ku9LG4UTUwUQcP1MiRJo65OyHh6Zn6huiEz78/Mv6NYDVSTWL26aMlY\ntmzQlUiS1F91Bn7eOcl7V06tnNFnl4kkaVx0eu+SsyJi5w73XRARfz61skaba2ZIksZBp/cuuRP4\nbkT8D3AxcA2wBrgPeDTwVOB5wCHAz4Ejel/q6PB+JpKkcdBRyMjMf4yI9wNvAI6iCBVVvwEuA47I\nzC/2tsTRVG3NMGRIkkZRpy0ZZOatwLuBd0fEo4EnAFtStHL8ODOzPyWOJlszJEmjrtMxGZ+MiG3K\n54cC92TmdZl5VWb+yIBRj2MzJEmjrNPZJQcBjyifnwPM6k8548WZJpKkUdZpd8mNwHsj4itAAK+J\niF+32jEzz+9VcePAsRmSpFHVacg4kuImaC8DEnhX+WezBAwZXXBshiRpVHU6u+QbwD4AEfEQ8AeZ\neXs/CxsntmZIkkZRnWXFZwN39LqQcebYDEnSKOp4CmtDZv6kH4WMu0ZrhiRJo6JOS4b6yOmskqRR\nYcgYEnaZSJJGjSFjiLg4lyRplHQdMiJiy4jYqvJ614j4m4g4oLeljR9bMyRJo6ROS8ZngEMBIuJR\nwNXA3wKfiYg39bC2sWRrhiRpVNQJGXOBr5XPDwZuA3alCB5H96iusWVrhiRpVNQJGVtR3Nod4ADg\nk5n5EHAVRdjQFNmaIUkaBXVCxo+AP4mIXYADgS+V27cHWt7PRN2xNUOSNArqhIwTgX8BbgKuzsxv\nltsPAL5dt5CIOCoiVkfEvRFxVUQ8ayP7z4qID0TEz8tjboyIP6r7+cPG1gxJ0kzXdcjIzE8ATwCe\nCVR/qV8OLK5TREQsAE4BlgB7AdcBl0bEdm323wy4rKzjlcDuwOHAz+p8/jCyNUOSNNN1vaw4QGbe\nCtwKEBHbAPsB38/MG2vWsRg4u3Gb+Ig4kuKOr4cB72ux/xuARwH7ZOaD5baba3720PLGaZKkmazO\nOhkXRcSby+dbAtcAFwHXR8SrapxvM2AeRUsIAJmZFC0V+7Y57OXAN4EzI+LWiPjfiHhbRIzU4mK2\nZkiSZrI6v5QnWDeF9U+BoGhVOBr4xxrn2w7YlGIqbNVtwI5tjpkDvJqi/pcAJ1Gs1XF8jc8fao3W\njOXLB1uHJEndqtNdMgu4q3z+R8B/Z+Y9EfE54OSeVVaEl2zz3iYUIeSIstXj2xHxeOAtwLsmO+ni\nxYuZNWvWetsWLlzIwoULp15xH8yeDRdcAIccAjvtBBMTg65IkjTKli5dytKlS9fbtnbt2lrnqhMy\nbgH2jYi7KELGIeX2RwP31TjfncCDwA5N27dnw9aNhjXA78qA0XADsGNEPCwzH2j3Yaeddhpz586t\nUebg7L33oCuQJI2LVl+8V65cybx587o+V53ukn8F/hP4KfBz4Kvl9gngf7s9WWbeD6wA9m9si4go\nX3+jzWH/AzypadvuwJrJAsZM5dgMSdJMVGcK65nAPhQzP55XrvYJsIp6YzIATgWOiIhDI+IpwFkU\nK4ueCxAR50fEeyr7fxDYNiL+LSKeHBEvA94GvL/m5w89182QJM00daewrqBofahu+1zdIjLzonJN\njBMpuk2uBQ7MzDvKXXYGHqjs/9Pyrq+nUayp8bPyeavpriOh2pqxapVTWiVJw69WyIiInYFXUCyG\n9fDqe5l5bJ1zli0kZ7Z5b78W264GnlPns2Yq182QJM0kXYeMiNgfuJiie+QpwHeAJ1LMBlnZy+K0\nPlszJEkzSZ2Bn+8F/iUzn04xm+RVwC7AlcB/9bA2teDYDEnSTFEnZOwBnF8+fwDYMjP/D3gH8NZe\nFabWnGkiSZop6oSMu4HNy+drgN0q77W8oZl6y1VAJUkzQZ2Bn1cBzwW+B3weOCUink5xN9Srelib\n2qiuAgqwYMFg65EkqZU6IeNYYOvy+ZLy+QLgh+V7mgaNYHHIIcWKoA4ClSQNm65DRmauqjy/Gziy\npxWpYy43LkkaZlO6NXpEbB0R21QfvSpMnXOmiSRpGHUdMiJidkR8LiLuBtYCvywfvyr/1DRxpokk\naZjVGZPxMYqFtw6juEtqu9uxaxpUZ5o4LkOSNEzqhIw9gXmZ+f1eF6PuOdNEkjSs6oSMb1Gs8GnI\nGBLONJEkDaM6IeONwFkR8XiK+5bcX30zM6/vRWHqTmOmiTdPkyQNizoh47EUq3yeU9mWFOM0Eti0\nB3WpS948TZI0bOqEjI8C3wYW4sDPoeKt4CVJw6ROyNgVeEVm/qjXxWhqqq0ZV14JExODrkiSNM7q\nLMZ1BcUMEw2hRmuGa2dIkgatTkvGZ4HTypui/S8bDvy8uBeFqZ5qa4bdJpKkQaoTMs4q/3xHi/cc\n+DkEJiYcBCpJGryuu0syc5NJHgaMIVEdBCpJ0iBM6QZpGn6OzZAkDYohY0Q1xmaArRmSpMEwZIyw\n6tgMWzMkSdPNkDHiqndplSRpOhkyRlz1Lq0XXjjoaiRJ46TrkBERX42IQyNiy34UpN5bsGBd0LDb\nRJI0Xeq0ZFwH/Atwa0R8OCL26XFN6oPGXVrtNpEkTZc662QcAzwOeD2wPbAsIr4XEW+JiB16XaB6\nw24TSdJ0qzUmIzMfyMxPZuYfAzsDHwdOAm6JiE9HxH69LFK9YbeJJGk6TWngZ0TsDZwAvAW4HXgv\ncCdwSUT8y9TLU6/ZbSJJmi51Bn5uHxF/GxHfAb4GPBY4BHhiZi7JzDcCBwBH9rZU9UK122TZskFX\nI0kaZXVukPZT4MfAR4FzM/OOFvtcD3xrKoWpfxqtGZIk9VOdkLF/Zn5tsh0y89fAC+uVpH6r3g7+\nyiuLlUElSeq1OmMyToiIRzVvjIhtIuKKHtSkadBYCdQlxyVJ/VInZMwHHt5i+xbA86dWjqaLN1CT\nJPVbx90lEfGHjafAUyNix8rbmwJ/BPysh7Wpz6o3ULPbRJLUa92MybgWyPLRqlvkXmBRL4rS9Kl2\nm6xaVbRwSJLUC910l8wGdqNoydi7fN14PB7YJjM/2vMK1Vd2m0iS+qXjlozM/En51Du3jhi7TSRJ\n/dBRyIiIVwBfyMz7y+dtZebFPalM08puE0lSr3XakvFpYEeKpcM/Pcl+STEIVDNMde2M5csNGZKk\nqesoZGTmJq2ea7RMTKxbcnynnew2kSRNjYFB62ksOb5mzWDrkCTNfJ2OyTi60xNm5un1y9GgVW+g\nBsXt4SVJqqPTMRmLO9wvAUPGDNcIFnabSJKmotMxGQ4DHDONbhNnm0iS6nJMhlqqLtK1fPlga5Ek\nzUydjsk4FXh7Zt5dPm8rM4/tSWUauOpsE3B8hiSpO52OydgL2KzyvJ2cWjkaNo7PkCTV1emYjBe2\neq7x4PgMSVIdUxqTERG7RMTOvSpGw8mbqEmS6ug6ZETEwyLipIhYC9wE/CQi1kbEuyJis40cPtl5\nj4qI1RFxb0RcFRHP6vC4QyLioYj4ZN3P1sZV722yevVga5EkzQx1WjLOAI4A/p5ifMZe5fM3UHON\njIhYAJwCLCnPdx1waURst5HjdgVOBpbV+Vx1ztkmkqRu1QkZfwb8ZWaenZnXl4+zgcPK9+pYDJyd\nmedn5o3AkcA95TlbiohNgI8B7wD8bj0NqrNNLrxw0NVIkoZdp7NLqn5L0U3S7Cbgd92erOximQe8\np7EtMzMiLgP2neTQJcDtmXlORDjnYZo420SS1Kk6LRnvB94eEZs3NpTP/6F8r1vbUdwe/ram7bdR\n3F5+AxHxXOD1wBtrfJ6mqDrbZJkdVZKkNjpdjKt5UOWLgJ9GxHXl6z2BhwOX97C2oMW6GxGxNfAf\nwOGZ+csefp461BifMX++01olSe112l2ytun1fze9nsrExjuBB4EdmrZvz4atGwC7AbsCn42IKLdt\nAhARvwN2z8y2YzQWL17MrFmz1tu2cOFCFi5cWK/6MTUxsS5oLF9uyJCkUbF06VKWLl263ra1a5tj\nQGcic/CLdEbEVcDVmXlM+TqAm4HTM/Pkpn0fDjyp6RTvBrYGjgZ+mJkPtPiMucCKFStWMHfu3D78\nFOPpwguL8RkXXOCy45I0qlauXMm8efMA5mXmyk6PqzPwsx9OBc6LiBXAcorZJlsB5wJExPnATzPz\n+Mz8HfC96sER8SuK8aI3TGvVWm8g6N5726IhSVqnVsiIiIOB1wBPoBiL8XuZ2XUzQWZeVK6JcSJF\nt8m1wIGZeUe5y87ABq0TGg477VT8ecklsGjRYGuRJA2POit+Hg2cQzFeYi+KlodfAHOAL9QtJDPP\nzMwnZuaWmblvZl5TeW+/zGy7ZkZmvj4zX1n3szU1ExNw+ulw9NGunyFJWqdOS8ZfA0dk5tKI+Evg\nfZm5KiJOBB7T0+o0YyxaBNtv7/oZkqR16qyT8QTgG+Xze4FHls//A3CKxhhz/QxJUlWdkHErsG35\n/GZgn/L5bIq1LTSmqvc3MWhIkuqEjCuAl5fPzwFOi4gvAxcCn+pVYZqZGutngHdslaRxV2dMxhGU\n4SQzPxARvwCeA1wMnN3D2jRDuVCXJAlqhIzMfAh4qPL6AuCCXhalma96x9bbb3dqqySNo7rrZDwa\neAOwB8X9RW4AzsnMu3pYm2a4BQuKgHH00cXME1cElaTxUmedjAlgNcUS3o+mmLZ6NLDaW66r2UEH\nFX8ecogDQSVp3NQZ+PkB4CJgdma+slwEaw5Fl8kHelmcZr7mGScOBJWk8VEnZDwJOCUzH2xsKJ+f\nyoY3LpPWm3GyfPlga5EkTZ86IWMlxViMZnsA102tHI2q6kDQM84YdDWSpOnQ0cDPiPjDysvTgX+L\niCcBV5Xb9gGOAo7rbXkaJdWBoOCME0kadZ3OLrmWYhZJdUXP97XY7+MUi3JJLTWChUFDkkZfpyHD\n5ZTUMwY7aGEaAAAQ0klEQVQNSRoPHYWMzPxJvwvReKkGDdfQkKTRVGfgJxGxW0ScERGXRcSXI+L0\niNit18VptLmGhiSNtjqLcR0IfA/YG7ge+A7wbOC7EfHi3panUeZdWyVptNVpyfgn4LTMfHZmHpuZ\nizPz2cC/Av/c2/I06prv2mrQkKTRUSdk7AF8pMX2jwJPnVo5GkcGDUkaTXVCxh3AM1psfwZw+9TK\n0bgyaEjS6KkTMj4MfCgi3hoRz4+I50XEccDZwId6W57GSXPQ8D4nkjSz1QkZJwEnAouAK4FlwJuB\ndwLv7lllGkuN5ccBLrlksLVIkqam08W4AIiIAHYBPpiZp0XEIwEy8zf9KE7jyeXHJWk0dBUyKJYV\n/xHw/4AfGi7UL64KKkkzX1chIzMfiogfAtsCP+xPSVLBoCFJM1u3LRlQ3Gn15Ih4U2Z+p9cFSVUG\nDUmaueqEjPOBrYDrIuJ3wL3VNzPzMb0oTGowaEjSzFQnZCymuO27NG2qQeMXv4B3vnOg5UiSOtB1\nyMjMc/tQh7RRixYVAeOEE4rXBg1JGm4dh4yI2AT4O+AVwMOBy4ETMvPeSQ+UeqgRLE44Abbd1q4T\nSRpm3bRkHE+x4NblFOMwjgF2AF7f+7Kk9t75ziJgOEZDkoZbNyHjdcBfZ+aHACLiRcDnIuINmflQ\nX6qT2mgeDHrQQcWt4yVJw6ObZcWfAHyh8SIzL6MYAPq4XhcldWLRIjj99CJozJnjTdUkadh0EzIe\nBtzXtO1+YLPelSN1Z9Gidfc68e6tkjRcuukuCeDciPhtZdsWwFkRcXdjQ2a+slfFSZ1YsAB22qkI\nGfPnF3dynZgYdFWSpG5aMs4DbgfWVh4fA37etE2adtW7t9qiIUnDoeOWjMx0FomG2oIFxZ+HHFIE\njQsuWLdNkjT9umnJkIbeggVFdwkUYeOMMwZbjySNszrLiktDbWICVq2CSy5xiqskDZItGRpJs2c7\nxVWSBs2QoZHmFFdJGhxDhkbeggUGDUkaBEOGxkJz0LjwwsHWI0njwJChsdFq5snq1YOtSZJGmSFD\nY6Ux88QBoZLUf4YMjZ3GzBPHaUhSfxkyNLYcpyFJ/eViXBpr1aXIDzkEbr/dhbskqVdsydDYW7Bg\nw3EatmpI0tQZMiQ2HKfh7BNJmjpDhlTRqlXDQaGSVM/QhIyIOCoiVkfEvRFxVUQ8a5J93xgRyyLi\nrvLx5cn2l7pRve8JOChUkuoaipAREQuAU4AlwF7AdcClEbFdm0PmAx8HXgDsA9wCfCkidup/tRoX\ndp9I0tQMRcgAFgNnZ+b5mXkjcCRwD3BYq50z8y8y86zMvD4zfwC8keJn2X/aKtZYcFCoJNU38JAR\nEZsB84DLG9syM4HLgH07PM0jgM2Au3peoMaeg0IlqZ6BhwxgO2BT4Lam7bcBO3Z4jn8GfkYRTKS+\ncFCoJHVnGEJGOwHkRneKOA54DfAnmfm7vlelseagUEnq3DCs+Hkn8CCwQ9P27dmwdWM9EfEW4O+B\n/TPzu5182OLFi5k1a9Z62xYuXMjChQs7LlhatAi23379lUL33LO4AZskzWRLly5l6dKl621bu3Zt\nrXNFMfxhsCLiKuDqzDymfB3AzcDpmXlym2P+DjgeOCAzv9XBZ8wFVqxYsYK5c+f2rniNtdWr4ZJL\niu4TKFo4XJZc0qhZuXIl8+bNA5iXmSs7PW5YuktOBY6IiEMj4inAWcBWwLkAEXF+RLynsXNE/D1w\nEsXsk5sjYofy8YjpL13jrNF9cuWVzkCRpGbD0F1CZl5UrolxIkW3ybXAgZl5R7nLzsADlUPeRDGb\n5BNNpzqhPIc0rSYmikdzF4qtGpLG2bC0ZJCZZ2bmEzNzy8zcNzOvqby3X2YeVnk9OzM3bfEwYGig\nFiwoWjXAVg1JGpqQIY2KiYliqmvzuhpOd5U0boaiu0QaNbNnF4+9915/YOiSJfC619mFImk82JIh\n9VF1YOiSJXDCCXahSBoftmRI06AxMHTbbYtWDdfWkDQODBnSNFq0qJhx0ry2hmFD0igyZEjTrNGF\nsueecN11jteQNLockyENyMRE+/Ea3uFV0iiwJUMasFbjNaCYArv33rZsSJq5bMmQhsSiRRuur2HL\nhqSZzJAhDZHZs4tVQw0bkkaBIUMaQtWwcfrpxTbDhqSZxpAhDbHqYl6tWjZcqlzSMDNkSDPAxETr\nbpT584v7otiyIWkYGTKkGaTVmI3q3V5t2ZA0TJzCKs1A1RuwLV9etGo0pr42VhDdZRenv0oaLEOG\nNINVw8Ytt6y/gii41oakwTJkSCOgETYmJtZfrry6sNdOO3l/FEnTy5AhjZjGCqJ77glr1tiVImlw\nDBnSiGq0WtiVImlQDBnSiOu0K8XWDUm9ZsiQxshkXSng2A1JvWXIkMZQc1eKYzck9YMhQxpjja4U\nWH/NjeaxG3anSKrDkCEJ2HDNDbtTJE2VIUPSeppbN6ozU5q7U8AWDkntGTIktdU8MwU2nAoL60KH\nLRySqgwZkjrSCBDVwNHoUmmEDls4JFUZMiR1rdpi0RgwevvtG7ZwNO4U66Jf0ngyZEiakuoYjuYW\njlaDRsFuFWlcGDIk9UxzC0erQaOwrltlzRpbOaRRZsiQ1BfVQaMHHVQEDmg9cLTayuFYDml0GDIk\n9V21S6XVTJVqKwesP4DUrhVp5jJkSJp21Zkq1VaO5tkqsH7gqB4rafgZMiQNVLWVA9bNVtlpp9Zd\nK0uWwH77rXtt6JCGlyFD0lBp17UCcMUVcMIJxaOh2tKxZg0sWDB9tUqanCFD0lCrtlRMTMAee6wb\nJNqqpeOGG9Zv6Wg+h6TpY8iQNKNUWyo6aemAorVj++2dwSJNN0OGpBmtuaWjuRWjVWsHrJs2u2aN\nt7KX+sWQIWmkNHeNNFo7GmGi1WqkDdX1OsAxHtJUGTIkjbzm4NFYjRSKIAHtg0erMR6tzilpQ4YM\nSWOnedosrB88GtqN8YANp9I2GD6kdQwZkkTr4NFqjAdMHj6aB5lWzyWNG0OGJE2iVThonkrb0G6Q\nKbQPH+0+QxoFhgxJqqHVgNDmQaYNk4UPmDyAOOtFM5khQ5J6qF3LR6vwARsPILDh/Vs6+UxpGBgy\nJGkatAsCkwWQVjeMa6VVCKmez9YQDYohQ5IGbLKWiOoN45p1GkJg/TVAWgUaMIyo9wwZkjTEWs16\nqWo19RY6W3yslXZdM83BxC4adcKQIUkz2MZCCGwYRKbSNdPQfPfbVi0j1fdsJRlPhgxJGnGdBBHY\neNfMZHe/7USrZds39llgq8lMZsiQJAGdh5Hmu99urCUDuuuyadZtq0m7mjX9hiZkRMRRwFuAHYHr\ngEWZ+a1J9n81cCLwROAHwHGZ+YVpKFVdWLp0KQsXLhx0GWPFaz79xvGad/tLu9XYkU5aMtq3miwF\nOr/m7ZaBrxtawODSiaEIGRGxADgFOAJYDiwGLo2IP8jMO1vsvy/wceCtwOeAPwM+HRF7Zeb3pq9y\nbcw4/s930Lzm089rvnGdtpI0a9dq8ra3LeW9713/mrcLBZMtAz8VnQ6S3dj2jb03k8ezDEXIoAgV\nZ2fm+QARcSTwMuAw4H0t9j8G+EJmnlq+XhIRBwBvBv56GuqVJE2TVi0G227beUtCu2Xgod4v/m4H\nyfZC83iWai29DDTt3vvBDzqvtWrgISMiNgPmAe9pbMvMjIjLgH3bHLYvRctH1aXAH/elSEnSjNZq\nGfipaDd1GHr7i7/RXVR3PMugDTxkANsBmwK3NW2/Ddi9zTE7ttl/x96WJknShup2/3Srubuo2XS2\nZBx++MbrbTYMIaOdALKH+28BcMMNN0ylJnVp7dq1rFy5ctBljBWv+fTzmk+/cbrmW2/d/r0nP7m7\n7XXf23LL3//u3KL90RsahpBxJ/AgsEPT9u3ZsLWi4dYu94diFgqvfe1ru69QUzJv3rxBlzB2vObT\nz2s+/bzmA/FE4Bud7jzwkJGZ90fECmB/4GKAiIjy9eltDvtmi/dfXG5v51Lgz4GbgPumVrUkSWNl\nC4qAcWk3B0VmNz0S/RERrwHOA/6KdVNYDwaekpl3RMT5wE8z8/hy/32BK4HjKKawLiyfz3UKqyRJ\nw2HgLRkAmXlRRGxHsbjWDsC1wIGZeUe5y87AA5X9vxkRC4F3l48fAn9swJAkaXgMRUuGJEkaPZsM\nugBJkjSaDBmSJKkvxiJkRMRREbE6Iu6NiKsi4lmDrmmURcTzI+LiiPhZRDwUEa8YdE2jLCLeFhHL\nI+LXEXFbRHwqIv5g0HWNsog4MiKui4i15eMbEfFHg65rnJT/7h+KiFM3vrfqiIgl5TWuProa+zjy\nIaNy87UlwF4Ud3i9tBxoqv54BMXg3aPobkE11fN84Azg2cCLgM2AL0XElgOtarTdQnGDxnnl4wrg\nMxGxx0CrGhPlF8XDKf5/rv76DsWEjB3Lx/O6OXjkB35GxFXA1Zl5TPk6KP4HcXpmtrr5mnooIh4C\n/iQzLx50LeOiDNC3AxOZ+fVB1zMuIuIXwFsy85xB1zLKImJrYAXwJuDtwLcz89jBVjWaImIJxczN\nuXXPMdItGZWbr13e2JZFqprs5mvSTPcoihakuwZdyDiIiE0i4hBgKyZfEFC98QHgs5l5xaALGRNP\nLru+fxwRH4uIXbo5eCjWyeijOjdfk2assqXuX4Gvu25Mf0XE0yhCxRbAb4A/zcwbB1vVaCvD3DOA\nZw66ljFxFfCXwPeBnYB3Assi4mmZeXcnJxj1kNFOtzdfk2aKM4GnAs8ddCFj4EZgT4qWo1cB50fE\nhEGjPyJiZ4oA/eLMvH/Q9YyDzKwuIf6diFgO/AR4DdBRt+Coh4w6N1+TZqSIeD/wUuD5mblm0PWM\nusx8AFhVvlwZEXsDx1CMFVDvzQMeC6woW+ygaKmeiIg3A5vnqA8yHLDMXBsRPwCe1OkxIz0mo0y7\njZuvAevdfK3ju8hJw64MGH8MvDAzbx50PWNqE2DzQRcxwi4Dnk7RXbJn+bgG+BiwpwGj/8pBt7sB\nHX+JGfWWDIBTgfPKO702br62FXDuIIsaZRHxCIqk2/i2MSci9gTuysxbBlfZaIqIMyluEvgK4O6I\naLTcrc1M7zjcBxHxbuALFDPVHklxh+f5wAGDrGuUlWMA1htnFBF3A7/IzBsGU9Voi4iTgc9SdJE8\nHjiB4j5iSzs9x8iHjA5uvqbeeybwFYpxL0mxTgkUd9o9bFBFjbAjKa7zV5u2vx44f9qrGQ87UFzb\nnYC1wPXAAc54mHa2XvTXzsDHgW2BO4CvA/tk5i86PcHIr5MhSZIGY6THZEiSpMExZEiSpL4wZEiS\npL4wZEiSpL4wZEiSpL4wZEiSpL4wZEiSpL4wZEiSpL4wZEiSpL4wZEiSpL4wZEiSpL4wZEgaiIjY\nLiLWRMRxlW37RsRvI+KFg6xNUm94gzRJAxMRLwE+DewLfB+4DvhUZv7dQAuT1BOGDEkDFRFnAC8G\nrgGeBjwrM+8fbFWSesGQIWmgImIL4DvAzsDczPzegEuS1COOyZA0aLsBj6P4/9HsAdciqYdsyZA0\nMBGxGbAc+DbFmIxjgadl5h0DLUxSTxgyJA1MRJwMvBL4Q+Ae4KvArzPz5YOsS1Jv2F0iaSAiYj5w\nNPDazLw7i288hwLPi4i/Gmx1knrBlgxJktQXtmRIkqS+MGRIkqS+MGRIkqS+MGRIkqS+MGRIkqS+\nMGRIkqS+MGRIkqS+MGRIkqS+MGRIkqS+MGRIkqS+MGRIkqS+MGRIkqS++P/cZzbfdo4nqQAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb48ddb1160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def exponential_dist(x, lam):\n",
    "    ''' \n",
    "    x: iterable of x values\n",
    "    lam: rate parameter value for the exponential distribution\n",
    "    returns: list of values of the PDF of an exponential distribution evaluated at every point in x\n",
    "    '''\n",
    "    return [lam*exp(-1*lam*point) for point in x]\n",
    "\n",
    "# lambda value\n",
    "l = 1\n",
    "_ = plt.plot(x, exponential_dist(x,l),'.', markersize = 1)\n",
    "_ = plt.xlabel(\"x\")\n",
    "_ = plt.ylabel(\"Probability mass f(x)\")\n",
    "_ = plt.title(\"Exponential distribution\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Long tails\n",
    "\n",
    "A distribution has a long tail if:\n",
    "\n",
    "$ \\lim_{x \\rightarrow \\infty} \\Pr[ X > x + t \\,|\\, X > x ] = 1 \\;\\; \\forall t > 0$\n",
    "\n",
    "Thinking about the simple case of a Power Law distribution:  \n",
    "CDF: $\\Pr[X > x] = cx^{-a}$  \n",
    "PDF: $p(x) = cax^{-(a+1)}$ \n",
    "\n",
    "We can solve the limit in the definition of long-tails:\n",
    "\n",
    "Definition: $ \\lim_{x \\rightarrow \\infty} \\Pr[ X > x + t \\,|\\, X > x ] $  \n",
    "\n",
    "Bayes law: $ \\lim_{x \\rightarrow \\infty} \\frac{\\Pr[ X > x + t \\, \\cap \\, X > x ]}{\\text{P}[ X > x]} $  \n",
    "\n",
    "Now $ X > x + t \\implies X > x $, so $\\Pr[ X > x + t \\, \\cap \\, X > x ] = \\text{P}[ X > x + t] $\n",
    "\n",
    "Rewritten: $ \\lim_{x \\rightarrow \\infty} \\frac{\\Pr[ X > x + t]}{\\text{P}[ X > x]} $\n",
    "\n",
    "Evaluate with the CDF: $ \\lim_{x \\rightarrow \\infty} \\frac{c(x+t)^{-a}}{cx^{-a}}  = \\lim_{x \\rightarrow \\infty} \\frac{(x+t)^{-a}}{x^{-a}} = 1$\n",
    "\n",
    "An intuitive explantion of this property from the wikipedia article: if the long-tailed quantity exceeds some high level, the probability approaches 1 that it will exceed any other higher level.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Subexponential distributions\n",
    "\n",
    "A subexponentail distribution satisfies what is called \"the catastrophe principle.\" To use a Twitter example: if Skippy, Josh, Brian, Jeff, and Kim Kardashian have 50 million Twitter followers between them, which scenario is more likely: \n",
    " 1. Each person has 1 million followers\n",
    " 2. One person has ~50 million followers (it's Skippy for sure)\n",
    "The answer is \"2\": one of the users has a huge number of followers, the others have a comparatively small number. This is called the \"catastrophe principle\"--one catastrophically large value increases the entire sum.\n",
    "\n",
    "The formal definition of a subexponential distribution is that it satisfies:\n",
    "$\\Pr[X_1+ \\cdots +X_n>x] \\sim \\Pr[\\max(X_1, \\ldots,X_n)>x] \\quad \\text{as } x \\to \\infty$\n",
    "\n",
    "In words: the probability that the sum of Twitter followers is 50 million is about the same as the probability that one of the users in the sum has 50 million followers.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pareto Distribution\n",
    "\n",
    "The Pareto distribution is (by definition) a Power Law Distribution.\n",
    "\n",
    "The PDF of the Pareto distribution is parameterized as: \n",
    "$\n",
    "f(x) = \\big \\{ \\begin{array}{ll} \\alpha \\, x_m^\\alpha \\, x^{-(\\alpha + 1)} & x\\geq x_m \\\\ 0 & x\\leq x_m \\\\ \\end{array} \n",
    "\\;, \\alpha \\geq 0 $\n",
    "\n",
    "The Pareto distribution is a heavy tailed distribution; it is also a long-tailed distribution, a power law distribution, and a subexponential distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAGHCAYAAAAdnkAlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xu4XGV5///3DUQQpVHhF0A5mKjFeFHRpETxEFQEPJeq\nJYn6Fc8FEfyGb7VoqwE86xeioCjtr5zKzw1YxSKtUgEFLcbYRKJoEAtBonIUG5GDHHL//lhrzGQy\ne2dm7TnP+3Vdc+3Za9bM3HvIxf7s57mfZ0VmIkmS1Gnb9LsASZI0mgwZkiSpKwwZkiSpKwwZkiSp\nKwwZkiSpKwwZkiSpKwwZkiSpKwwZkiSpKwwZkiSpKwwZknouIg6MiI0RsbAH73VCRGxsOLYxIk7t\n9nuX7/Wm8v326sX7SYPEkCF1QUQcUf5iqd3ui4ifRcRpETGrTzU9MiKWdfoXe0Ts3fCzPhARd0TE\nf0bERyJiz0me2vY1DSLifRHxF20+LYGNWz1rmqaoLanws0qjILx2idR5EXEEcCbwAeAmYAfgecAb\ny+/3zcz7e1zTzsAdwAmZeVIHX3dvYB3wReDfKf54eSywP/Dq8rS3ZuYFDc97RGY+0OZ73Q18KTPf\n0sZztgG2q3+vcmTjs5l5bDvvX6W2iAhgRrs/qzQKtut3AdKI+0Zmri7vnxkRdwFLgb8ALpj8aVsX\nEdsDD2TrfynEdN6vBasz84ubvWHE3wHfBM6OiJ9m5o9rj3X7l25E7JiZ92bmRqBvv+DL/z4GDI0l\np0uk3rqC4pf9bICIeGxE/N+I+FFE3B0RGyLi3yPi6fVPquthWBQRH46I9cA9wE7l4zMj4tMRcXNE\n3B8RP4+I95Z/RddGG26nGLY/oW5q44N17/GiiPhORPw+In4bEV+NiKdO54fNzPXAm4Dtgfc2+XkW\n1h17ckR8OSJuKaeX1kfERETUfsaNwI7Am+rqP7N8rPYzzY2IL5Zh7jv1jzWrLyJeFxHXle/3XxHx\n/IbHz46IdU2et9lrbqW2pj0ZEfHOiLi2/O/1q4j4bETMbDjn2+W/jbkR8a2IuCcifhkR79nqhy8N\nAEcypN56cvn1N+XXOcCrgC9RTDnsCvw18O2IeFpm3trw/A8AfwD+L8Uv7gci4pHAVcDjgc8D64Hn\nAB8DdgOOo5gmORL4AvCV8gbwI4CIeDHFVMcNwDLgkcCxwHcjYl5m3lz1B87MFRFxA3Bw40O1OxEx\nA/gPYAZwKnAr8ATgFcBjgLuBNwD/BHwf+IfyqTc0vNaXgOuB97Fp5GaynogXAIvK9/sD8E7g6xGx\nIDN/upXnNh7fWm2bvUZEnAB8sPyZTwf2Kd//zyPiuZn5cN1zHwd8neK/2fnAa4GPR8SPMvPSJrVJ\ngyMzvXnz1uEbcATwMPBCYGeKX5iLKH7Z/x7YvTxvRpPn7gXcB/xd3bEDKZoXfw48ouH8vwd+B8xp\nOP5RimH6J5Tf71y+xgebvOcPgVuAmXXH/gx4CDhrKz/r3uXrHjfFOReVn8ej636eh4GF5ff7la/x\nl1t5r7uBM5scX1Y+/7xJHnu44djG8v2fUXdsT+Be4F/qjp0F3Njia05WW+3fwl7l97sA9wP/3nDe\nO8vzjqg79q3y2Ovqjs0o/1td2O9/5968be3mdInUPQFcThEs1lM0Rv4OOCwzbwHIzAf/eHLENhHx\nOIpfdD8D5jV5zbNzy16G11JMDWyIiJ1rt/K9twOmXE0SEbtR/JI/KzM31I5n0T/xTeBlrf/Ik/p9\n+XWnSR6vve9LypGZKpJipKZVV2fmNX98cjG186/AIbVppi55MUVQ+HTD8X+kCCovbzh+T9b1upT/\nZr5PMQomDTRDhtQ9CRxF8UvlBcDTMvNJmXlZ7YQoLI2I6ymG7O+k6J34M2Dmli/JTU2OPQV4CUWY\nqb99s6xha0tm9y6/Xt/ksbXALtP4xV/z6PLr3c0ezMybgJOBtwF3RsQ3yp6FP2nzfbbon5jCfzc5\ndj3wKIrRhm5p+nmX4eHGusdr1jd5jd9SrOCRBpo9GVJ3/SA3rS5p5u+Akyjm8/8euItiKP8zNP8j\n4L4mx7ahCBSfoPkKkmbhoV63V50A7Avcnpm/n+yEzHxPRJxNsfLmEIpeieMj4tmZ+esW36fZ59OO\nxs9ispU723bwPbbm4UmO9+K/mzQthgypv14DXJGZb68/GBGPoRiNaMUNFL0O39rKeZP9wryp/LpP\nk8eeCtyZmZV/eUfEAcCTgHO3dm5m/gT4CfDRiHg2cDVFw2ptFUwnN/Z5SpNjf0oxXXVn+f1vKRpP\nGz2xybFWa7up/LpP3f1a8+tsisAojQSnS6T+epiGv0gj4q8oGkVbdSFwQEQc0vhAubS19lf3veXX\nzX5pZrGC5RrgiPrpiYjYl2JE4d/aqKXx/fcGzmbTipjJztuprs6an1CM6mxfd+yexvqn4YCI+GPf\nSxQ7k74KuDQza4HhBmBm+VnUztsdOKzJ67Va22XAgxSrd+q9DfgT4JKWfwJpwDmSIXVPK8PZlwAf\nKPdUuJqiF+P1bFr+2IpPUfxyvKScblhF0VfwdIodN58I3JWZ90fET4FFZQ/Ib4Fry9GD91AsYV0R\nEf9EsefDu8pzTmyxjvkR8XqKP14eQ7Hj52sogsIbMvPahvPrP58XAZ+NiNoS1O0odkd9CPhy3Xmr\ngBdHxFLg18C6zFzZYn2NrqVYsnoaxSqcoyj3Eak7Z4JiGuqrUVzr5FEUIyvNGnNbqi0z74yIjwEf\njIhvABdTjBgdBawE/r+KP480cAwZUve0Mnz+UYpf6K8DDqf4RfUy4ONNnt/09TLzvnJTq/cDfwX8\nL4pVLNdTTDNsqDv9rcBpwHLgERQB4ieZeXlEvKT8/kSKv7S/DRyfmb9o4edIYHF5e6h8/58DpwBn\nZOYvJ3lOzRrgGxT7YjyBYtRlDfCShl/UxwFnAB+i2MvjHIpfzK3U1/j9t4EVFKFiT4qRkzfWh6HM\n/G1EHFb+HJ+gaCw9nmJapTFktFxbZp4YEbdTBLlTKHpxvkCxbLmxB2Oyf0deE0IDz2uXSJKkrhiI\nnoyIeH5EXFxurbsxIl7VwnNeEBGryi15r4/iglSSJGlADETIoJjnvAY4mhaGACPiiRRz2ZdTbCL0\nGeD/jYjGbYslSVKfDNx0SXmhocMy8+IpzvkE8NLMfHrdsQmKLZE7sTuhJEmapkEZyWjXsymWgdW7\nFDigD7VIkqQmhjVk7Abc1nDsNuBPImL7JudLkqQeG6UlrPWXdd7yweKCUYdS7LB3f49qkiRpFOxA\nsefOpZn5m1afNKwh41Zg14Zjs4DfNblCZc2huMmNJEnT8XqKK0q3ZFhDxveAlzYcO6Q8PpmbAM47\n7zzmzp3bpbLU6KijlvL5zy/vdxljZenSpSxf7mfeS37mvedn3ltr167lDW94AzS/EvSkBiJkRMSj\ngCezacpjTkTsR7EV8vpyC97HZ2ZtL4wvAO8qV5mcCRwEvJZip8TJ3A8wd+5c5s1r3KhP3bLTTjP9\nvHts5kw/817zM+89P/O+aavdYFAaP/8c+CHFlsoJnAysZtM1E3aj2PYXgMy8CXg58GKK/TWWAm/N\nzMYVJ+qjdevg8suLr5Kk8TMQIxmZeSVTBJ7MfPMkz5nfzbo0PbNnw0EHFV8lSeNnUEYyNKJ23NGR\nDEkaV4YMddVBBy1hzhyDRi8tWbKk3yWMHT/z3vMzHw4Dt614t0TEPGDVqlWrbBbqsXXrnDKRpGG2\nevVq5s+fDzA/M1e3+jxHMiRJUlcYMtRV69bhdIkkjSlDhrpq9my48UanSyRpHBkyJElSVxgy1FVO\nl0jS+DJkqKucLpGk8WXIUNfNnu1IhiSNI0OGus4pE0kaT4YMdZ1TJpI0ngwZkiSpKwwZ6jqnSyRp\nPBky1HVOl0jSeDJkqCdcYSJJ48eQoZ5wykSSxo8hQz0xezZceaVTJpI0TgwZ6ol16+DAAx3JkKRx\nYshQT9j8KUnjx5ChnrH5U5LGiyFDPWPzpySNF0OGesYpE0kaL4YMSZLUFYYM9YzTJZI0XgwZ6hmn\nSyRpvBgy1FOuMJGk8WHIUE85ZSJJ48OQoZ5ye3FJGh+GDPWU24tL0vgwZKinbP6UpPFhyFDP2fwp\nSePBkKGes/lTksaDIUM9Z/OnJI0HQ4Z6zuZPSRoPhgz1nM2fkjQeDBnqC5s/JWn0GTLUFzZ/StLo\nM2SoL2z+lKTRZ8hQX9j8KUmjz5ChvrD5U5JGnyFDfWPzpySNNkOG+sbmT0kabYYM9Y3Nn5I02gwZ\n6hubPyVptBky1DeOZEjSaDNkqG8cyZCk0WbIUN+4jFWSRpshQ33lMlZJGl2GDPWVy1glaXQZMtRX\nNn9K0ugamJAREUdHxLqIuC8iVkTE/ls5/39HxHURcW9E3BwRp0TE9r2qV51h86ckja6BCBkRsQg4\nGVgGPBNYA1waEbtMcv7rgI+V5z8VeAuwCPhITwpWx9SaPyVJo2cgQgawFDgjM8/NzOuAI4F7KcJD\nMwcA383MCzLz5sy8DJgAFvSmXHWafRmSNHr6HjIiYgYwH7i8diwzE7iMIkw0czUwvzalEhFzgJcB\n/9bdatUNLmWVpNG0Xb8LAHYBtgVuazh+G7BPsydk5kQ5lfLdiIjy+V/IzE90tVJJktSyvo9kTCGA\nbPpAxAuA91NMqzwTeDXwioj4+55Vp45xGaskjaZBGMm4E3gY2LXh+Cy2HN2oOQk4NzPPKr//SUQ8\nGjgD+PBUb7Z06VJmzpy52bElS5awZMmSdutWh9j8KUmDY2JigomJic2ObdiwodJr9T1kZOaDEbEK\nOAi4GKCcAjkIOHWSp+0IbGw4trF8apQ9HU0tX76cefPmTb9wddycOfZmSFK/NfvDe/Xq1cyfP7/t\n1+p7yCidApxTho2VFKtNdgTOBoiIc4FfZub7y/O/BiyNiGuA7wNPoRjd+NepAoYGl5tySdLoGYiQ\nkZkXlo2cJ1FMm1wDHJqZd5Sn7AE8VPeUD1GMXHwIeAJwB8UoiD0ZQ6q2KZcjGZI0OgYiZABk5unA\n6ZM89qKG72sB40M9KE09YF+GJI2eQV5dojHkKhNJGh2GDA0M+zIkabQYMjQwvFiaJI0WQ4YGhn0Z\nkjRaDBkaOPZlSNJoMGRooHixNEkaHYYMSZLUFYYMDRQvliZJo8OQoYFi86ckjQ5DhgaSoxmSNPwM\nGRo4bsolSaPBkKGB46ZckjQaDBkaOI5kSNJoMGRo4DiSIUmjwZChgeMKE0kaDYYMDSxXmEjScDNk\nSJKkrjBkaCA5ZSJJw6/tkBERR0xyfLuI+Nj0S5I2ccpEkoZXlZGM0yLiSxHx2NqBiHgqsBJY3LHK\nNPZcyipJw61KyHgm8ATgxxFxcEQcDawG1gLP6GRxGm8uZZWk4bZdu0/IzBsi4nnAp4FvAA8DR2Tm\nRKeL03izL0OShlvVxs9XUEyNXA38D/C2iHh8x6qS6tiXIUnDqUrj5xnAhcAngIXA04E/UEyfHN7Z\n8jTu7MuQpOFVZSTjucCzMvPkLNyamS8DPgic2dnyNO7sy5Ck4dV2TwYwPzP/0HgwMz8XEZd1oCbp\nj+zLkKTh1fZIRrOAUffYz6ZXjtScfRmSNHzc8VMDz74MSRpOhgwNvFpfxlVX9bsSSVI7DBkaeLWR\nDBtAJWm4VFnC+saI2L7J8UdExBs7U5a0uYULnTKRpGFTZSTjLGBmk+M7lY9JHedSVkkaPlVCRgDZ\n5PgewIbplSM151JWSRo+Le+TERE/pAgXCVweEQ/VPbwtMJviWiZS18yZU4QNp00kafC1sxnXV8uv\nzwAuBX5f99gDwE3AlztTlrQll7JK0nBpOWRk5okAEXETcEFm3t+toqRman0ZV15ZNIJKkgZblR0/\nzzFgqB9cyipJw6Xta5dExLbAUuBwYC/gEfWPZ+bjOlOatCWXskrS8KiyumQZcBxwAcVS1lOArwAb\ngRM6VpnUhEtZJWl4VAkZrwfenpknAw8BE5n5NuAk4NmdLE5qZPOnJA2PKiFjN+DH5f3fs2ljrkuA\nl3eiKGkyXsdEkoZHlZDxS2D38v4NwCHl/f2BSS8DL3WCzZ+SNDyqhIyLgIPK+6cBH4qInwPnAmd2\nqjBpMjZ/StJwaHt1SWYeX3f/goj4BfAc4OeZ+bVOFic1434ZkjQc2g4ZjTJzBbCiA7VILamfMnGL\ncUkaXFWmS6S+c8pEkgafIUNDyf0yJGnwGTI0lLz0uyQNPkOGhtqcOY5mSNKgajtkRMSeEbFH3fcL\nIuLTEfGOzpYmSZKGWZWRjC8CLwSIiN2AbwILgI9ExAc7WJs0JadMJGmwVQkZ+wIry/uHA9dm5nMo\nrmnypqqFRMTREbEuIu6LiBURsf9Wzp8ZEZ+LiF+Xz7kuIl5S9f01vJwykaTBVGWfjBls2j78xcDF\n5f3r2LTdeFsiYhFwMvAOigCzFLg0Iv40M+9scv4M4DLgVuDVwK+BvYH/qfL+kiSp86qMZPwEODIi\nng8cDHyjPP544DcV61gKnJGZ52bmdcCRwL3AWyY5/63AY4DDMnNFZt6cmd/JzB9Pcr5GlFMmkjS4\nqoSMvwX+Gvg2xWXe15THX8WmaZSWlaMS84HLa8cyMylGKg6Y5GmvBL4HnB4Rt0bEjyPifRHhapkx\n5ZSJJA2eKtcu+XZE7AL8SWb+tu6hf6AYfWjXLsC2wG0Nx28D9pnkOXOAFwHnAS8FngKcXr7OhyvU\nIEmSOqztkBERjwSiFjAiYm/gL4G1mXlpB2sLICd5bBuKEPKOctTjhxHxBOBv2ErIWLp0KTNnztzs\n2JIlS1iyZMn0K1ZfOGUiSZ0zMTHBxMTEZsc2bNhQ6bWqNH7+K/AV4AsR8Rjg+8CDwC4RcVxmfr7N\n17sTeBjYteH4LLYc3ai5BXigDBg1a4HdImK7zHxosjdbvnw58+bNa7NEDYM5c7xgmiRNV7M/vFev\nXs38+fPbfq0qPQzzgO+U919LEQT2Bt4IHNvui2Xmg8Aq4KDasYiI8vurJ3nafwJPbji2D3DLVAFD\no2/9+n5XIEmqqRIydgTuLu8fAnwlMzdSXO5974p1nAK8IyLeGBFPBb5Qvs/ZABFxbkR8tO78zwM7\nR8RnIuIpEfFy4H3AZyu+v4Zc/eXfbQCVpMFQZbrkv4HDIuIi4FBgeXl8FvC7KkVk5oVlM+lJFNMm\n1wCHZuYd5Sl7AA/Vnf/LiDikfO81wK/K+5+s8v4aDV7+XZIGS5WQcRLF1uLLgcsz83vl8UOAH1Yt\nJDNPp1gh0uyxFzU59n3gOVXfT6Ondvl3+zIkaTBUWcL6LxHxXYrdPdfUPXQ5cFGnCpMkScOtykgG\nmXkrxZbe9cfa3ohL6iSXskrSYKkUMsqLl/0VsBfwiPrHMvPVHahLqsylrJI0GNpeXRIRiymWkM6l\n2IRrBvA0ih04q+3WIXWYS1klqf+qLGF9P7A0M18JPAC8myJwXAjc3MHapLa5lFWSBkeVkPEk4N/K\n+w8Ajyp33lxOcal2qa9cyipJg6FKyLgL2Km8/ytg3/L+Yyg20JL6qraU9aqr+l2JJI23KiHjO8DB\n5f0vAZ+JiH8EJqi7XLvUL06ZSNJgqLK65F3ADuX9j1BcHO05wJfxMusaEHvuWXxdv95pE0nqlyqb\ncd1Vd38j8PGOViR1QP1ohstZJak/Ku2TARARsyiuV7LZlEtm/mi6RUmdYAOoJPVX2yEjIuYD51As\nW42GhxPYtgN1SdNWawC98soicEiSeqvKSMaZwPXAW4HbKIKFNHCcMpGk/qoSMuYAr8nM/+50MVKn\n2QAqSf1TZQnr5cB+nS5E6gaXs0pS/1QZyXgbcE5E7AtcS7GE9Y8y8+JOFCZ1ig2gktQfVULGAcBz\ngZc2eczGTw0cG0AlqT+qTJecBpwH7J6Z2zTcDBgaOE6ZSFJ/VAkZOwPLM/O2ThcjdUutAVSS1DtV\nQsZXgBd2uhCpF9av73cFkjQ+qvRkXA98LCKeB/yYLRs/T+1EYVInuWeGJPVe1dUlvwcOLG/1EjBk\naCC5ykSSeqvKBdL8X7SGkqtMJKm3qvRkSEPJVSaS1FuGDI2V+m3GJUndZcjQWHE0Q5J6x5ChseNo\nhiT1hiFDY8fRDEnqjbZDRkTcFBEfjIi9ulGQ1AuOZkhS91UZyfgM8Grgxoj4ZkQsjojtO1yX1FWO\nZkhS97UdMjJzeWY+A1gArKW4YNotEfHZiJjX6QKlbnE0Q5K6q3JPRmauzsxjgccDJ1HsBPqDiFgT\nEW+JiOhUkVI3OJohSd1VZVtxACJiBvCXwJuBg4EVwD8BewAfBV4MvK4DNUpd49VZJal72g4Z5ZTI\nm4ElwMPAPwNLM/O6unMuAn7QqSKlblu/3muaSFKnVRnJ+AHwTeAo4KuZ+WCTc9YB50+nMKkXvDqr\nJHVPlZAxJzN/MdUJmXkPxWiHNPDqG0ANGZLUOVUaP78VETs3HoyIx0TEjR2oSeopG0AlqTuqhIwn\nAts2Ob498IRpVSP1ictZJanzWp4uiYhX1X17aERsqPt+W+Ag4KYO1SX1lL0ZktR57fRkfLX8msA5\nDY89SBEw/k8HapL6wt4MSeqslkNGZm4DEBHrgP0z886uVSX1gaMZktRZVbYVn23A0KiyN0OSOqel\nkYyIOLbVF8zMU6uXI/WXoxmS1DmtTpcsbfG8BAwZGmr2ZkhSZ7QUMjLT/9VqbDiaIUmdUfkqrNIo\nszdDkqav1Z6MU4APZOY95f1JZeZxHalM6iNHMyRp+lrtyXgmMKPu/mRyeuVIg2flSkOGJFXRak/G\nC5vdl0bZwoVw/vmweDEsWGDQkKR2TasnIyL2jIg9OlWMNGgWLCi+2pshSe1rO2RExHYR8aHy2iU3\nAb+IiA0R8eGImLGVp0tDxSu0SlJ1VUYyTgPeAbyXoj/jmeX9tzKNPTIi4uiIWBcR90XEiojYv8Xn\nLY6IjRHxlarvLU3FlSaSVE2VkPE64E2ZeUZm/qi8nQG8pXysbRGxCDgZWEYRWtYAl0bELlt53t7A\np4Crqryv1Ir60Yyr/JcmSS2rEjL+QPNLut8EPFCxjqXAGZl5bmZeBxwJ3EsRXJqKiG2A84APAg5k\nq6tqoxlOm0hS66qEjM8CH4iI7WsHyvt/Vz7WlrKPYz5wee1YZiZwGXDAFE9dBtyemWe1+55Su2qj\nGeC0iSS1qtXNuBr7HV4M/DIi1pTf7wc8grqg0IZdgG2B2xqO3wbsM0k9zwXeXL6v1BMLF7pBlyS1\no9XNuDY0fP/lhu+78bdd0GRzr4h4NPDPwNsz87ddeF9pUrVpEzfokqSta3Uzrjd3sYY7gYeBXRuO\nz2LL0Q2AJwF7A1+LiCiPbQMQEQ8A+2TmpLPmS5cuZebMmZsdW7JkCUuWLKlWvcbK7NmbNugCWLSo\nv/VIUqdNTEwwMTGx2bENGxrHGloTRftDf0XECuD7mfnu8vsAbgZOzcxPNZz7CODJDS/xEeDRwLHA\nzzPzoSbvMQ9YtWrVKubNm9eFn0Lj5IILiqDhtImkcbB69Wrmz58PMD8zV7f6vFanSzYTEa8FDgf2\noujF+KPMrPIb/BTgnIhYBaykWG2yI3B2+X7nAr/MzPdn5gPATxvq+Z/irXNthfeW2la/E6ghQ5Ka\nq7Lj57HAWRRTGc+kCAW/AeYAX69SRGZeCPwf4CTgh8DTgUMz847ylD2A3aq8ttRNLmmVpMlVWcL6\nTuAdmXkMxb4Yn8zMgyl2+5w55TOnkJmnZ+YTM/ORmXlAZv5X3WMvysxJ98zIzDdn5qurvrfULpe0\nStLWVQkZewFXl/fvA3Yq7/8zYPekxkb9klZ3ApWkLVUJGbcCO5f3bwaeXd6fTbHsVBob7gQqSZOr\nEjKuAF5Z3j8LWB4R3wQuAC7qVGHSMKifNlm5sr+1SNKgqbK65B2U4SQzPxcRvwGeA1wMnNHB2qSh\nsHChe2dIUjNth4zM3AhsrPv+fOD8ThYlDZtasFi8uFje6rJWSaq+T8ZjgbcCcym2/l4LnJWZd3Ww\nNmmo1PbOcMtxSSpU2SdjIcWl1Y8FHgs8rry/rnxMGkv1W45fcEG/q5Gk/qsykvE54ELgqMx8GCAi\ntgVOLx/7s86VJw0Xp00kaZMqq0ueDJxcCxgA5f1T2PKaItLY2X334qurTSSNuyohYzVFL0ajucCa\n6ZUjDb/61SZu0iVpnLU0XRIRT6/79lTgMxHxZGBFeezZwNHA8Z0tTxpOtSbQAw8s9tFYaLeSpDHU\nak/GNRSrSOp39Pxkk/O+SLEplzTWapt0HXhgcfOS8JLGUashw/89Sm2qv7aJy1oljaOWQkZm/qLb\nhUijyN1AJY2zqptxPQn432y+GddnMvOGDtYmjYT6Za27725/hqTxUWUzrkOBnwILgB8B1wLPAn4S\nEQd3tjxpNNQ3gnq1VknjosoS1o8DyzPzWZl5XGYuzcxnAZ8GPtHZ8qTRUNsNFOCSS/pbiyT1SpXp\nkrnA4U2On0kxhSKpiUWL4Pbb4dhjYdYs+zMkjb4qIeMO4BnAzxuOPwO4fdoVSSPsmGOKgGF/hqRx\nUCVk/CPwDxExB7iaovHzecDfAid3sDZpJLlRl6RxUaUn40PAScAxwJXAVcC7gBOAj3SsMmlE1Tbq\ngiJouPW4pFHVVsiIiAD2BD6fmXsAM4GZmblHZn4mM7MbRUqjprZRFxg0JI2udkcyAvhviqBBZt6d\nmXd3vCppDDQGDZe2Sho1bYWMzNxI0fC5c3fKkcZLfdDw0vCSRk2VnozjgU9FxL6dLkYaR/Vbj1/g\n5QUljZAqq0vOBXYE1kTEA8B99Q9m5uM6UZg0TmorTlzaKmmUVAkZSymWrUrqkMZLw7u0VdIoaDtk\nZObZXahDGnv1l4Y3aEgaBS33ZETENhHxtxHxnxHxg4j4eEQ8spvFSePGpa2SRkk7jZ/vp9hs6/fA\nr4B3A6d3oyhpnDUGDZtBJQ2rdkLGEcA7M/PQzDwMeCXwuoioskJF0hTqg8bixXDaaf2tR5KqaCcg\n7AV8vfZNZl5G0QD6+E4XJakIGjfeCKeeWly51RENScOmnZCxHXB/w7EHgRmdK0dSvdmz4RWvKO4v\nXmyPhqTRQyI+AAAOMElEQVTh0s7qkgDOjog/1B3bAfhCRNxTO5CZr+5UcZJc3ippeLUzknEOcDuw\noe52HvDrhmOSOsxmUEnDqOWRjMx8czcLkTS1+n00Fi+G22+HY47pd1WSNLkqO35K6pNaM+gllxTN\noGDQkDS4DBnSkJk9e1OwMGhIGmTucSENqdqqk2OPdR8NSYPJkQxpSM2eveXUySteURyXpEHgSIY0\nxGpTJ7UNu+bMcS8NSYPDkCGNgGOOgfPPL+67xFXSoDBkSCNi0aLNr3dywgl9LUeS7MmQRklties5\n58CJJxbHDBuS+sWQIY2Y2bM3BYsTT4Sdd7YhVFJ/OF0ijagTTti8IdQ+DUm9ZsiQRlh9Q+jixcV+\nGuvW9bcmSePDkCGNuEWLij4NRzUk9ZohQxoDtf00Gkc1JKmbbPyUxsiiRbD77sVeGu4SKqnbHMmQ\nxkxtmavTJ5K6bWBCRkQcHRHrIuK+iFgREftPce7bIuKqiLirvH1zqvMlbW6y6RO3JJfUSQMRMiJi\nEXAysAx4JrAGuDQidpnkKQcCXwReADwbWA/8R0Ts3v1qpdHR2BR64IFu3iWpcwalJ2MpcEZmngsQ\nEUcCLwfeAnyy8eTM/F/130fE24DXAAcB53W9WmmE1EY19tsPrrhi006hRxxhr4ak6el7yIiIGcB8\n4KO1Y5mZEXEZcECLL/MoYAZwV+crlMbDwoXFDYqgceKJxXTKggWGDUnVDMJ0yS7AtsBtDcdvA3Zr\n8TU+AfwKuKyDdUlj6YQTNu/VsDFUUlWDEDImE0Bu9aSI44HDgcMy84GuVyWNgVqvRn3YOOEEdwuV\n1J6+T5cAdwIPA7s2HJ/FlqMbm4mIvwHeCxyUmT9p5c2WLl3KzJkzNzu2ZMkSlixZ0nLB0jiYPbu4\nLViw6aquTqFIo29iYoKJiYnNjm3YsKHSa0XmVgcLui4iVgDfz8x3l98HcDNwamZ+apLnvAd4P3BI\nZv6ghfeYB6xatWoV8+bN61zx0pg47bRNG3hBETYWLepfPZJ6Z/Xq1cyfPx9gfmaubvV5gzJdcgrw\njoh4Y0Q8FfgCsCNwNkBEnBsRf2wMjYj3Ah+iWH1yc0TsWt4e1fvSpfFwzDFbTqG4t4akqQzCdAmZ\neWG5J8ZJFNMm1wCHZuYd5Sl7AA/VPeUoitUk/9LwUieWryGpC+qnUC65ZNPIxrJlLnmVtKWBmC7p\nBadLpM676qrN99Y4//zi2ii1pbCSRkPV6ZKBGMmQNJxqe2vsvHMxqrF4cXH81FOLzb0MG9J4M2RI\nmrZjjimu5rp+PaxZs2kaxbAhjTdDhqSOqPVrLFxYBAvDhiRDhqSOq02jzJpVTKEYNqTxZMiQ1DWL\nFhUrUSabRtlzT1ekSKPMkCGpq6aaRgF3EJVGmSFDUs/UplHqw0ZtRYrLX6XRY8iQ1HP1YeOWW4qg\nUb/8ddYstyyXRoEhQ1Lf1EYtmvVtrF0LL3rR5udJGi6GDEl919i3UdtFtLaTqI2i0nAyZEgaKLWp\nlNooRmOjqMtgpeFhyJA0kGohorF3o3EZbP25kgaLIUPSwKvv3Vi5Em6/vfnoRv25kvrPkCFpaNR6\nN2BTqJhsOsX+Dan/DBmShtLWplNg094b9edL6h1DhqSh17gUtnHvDYBly1wSK/WaIUPSyKifTqkF\nDph8SSw4rSJ1kyFD0kiqDxwLF8LcucXUSWMPB2yaVjFwSJ1lyJA0FmrblNd6OGDL66fA5n0chg5p\negwZksZOY9MoNO/jAJtHpekwZEgaa/XBob6Po9koR30vR+NzJW3JkCFJpcY+jvpA0ayXw9AhTc2Q\nIUmTqA8NzaZWJgsdt9xSTLEYOjTuDBmS1KLJplagtZEOG0k1bgwZklRB/dQKbDnS0Xh9lZr6RtLa\n86RRZciQpA5pDAzNejoaV680TrE0ex1pWBkyJKlLJuvpqGk2xQJbTrPccsumfT6kYWLIkKQeaRyh\naJximWxHUoC1azdde6XGHg8NOkOGJPXRVMGjpvHaK/UaRz3A8KHBYciQpAHTLHg0jmJMtoy2pln4\nsNdDvWbIkKQh0CwgLFgAK1duvlplqvDRrMkUHPlQ9xgyJGlINS6jrWkWPibr9aipLa01gKiTDBmS\nNGKahY9mTaa1+80uDFfv/POLr/XhAwwg2jpDhiSNicl6Mup3L+1UAAFDiAwZkjT2Jpt2gU0BpHEa\npZUAAs0bUGsMIaPPkCFJmlSVAAJbX/1S07jNeu257n46GgwZkqRKpgog0LwBtWaybdYbNdv9tNnr\nGUYGkyFDktQVU4WQZpuO1Wxt99NmTj0VZs1qHkDqX9Mw0luGDElSX2ztF36zINJsJKPdMNLKyAjY\nM9IJhgxJ0sBqZeShFkamCgy33AK33956GKlp1jNSez2X9G6dIUOSNPRanQZpZWQEWu8ZaTTVkt5x\n7CcxZEiSxkarv9Cn6hmB6kt6m2m1n2SyOgeZIUOSpCba/QU+1ZJemH4/STPLlm158bzJ3qvZY90O\nKYYMSZI6YGtLeptptZ+k2WNXXAEnnljcqmq1Efb666u9viFDkqQ+qjqasHAhzJ3bXv9H/WPTHUVp\nhSFDkqQhtWhR9ee2ukQYipGMt7+9/fcwZEiSNKZaHUV59KOrvf421Z4mSZI0NUOGJEnqCkOGJEnq\nCkOGJEnqioEJGRFxdESsi4j7ImJFROy/lfP/KiLWlueviYiX9qpWtW5iYqLfJYwdP/Pe8zPvPT/z\n4TAQISMiFgEnA8uAZwJrgEsjYpdJzj8A+CLwj8AzgK8CX42Ip/WmYrXK/xH0np957/mZ956f+XAY\niJABLAXOyMxzM/M64EjgXuAtk5z/buDrmXlKZv4sM5cBq4F39aZcSZK0NX0PGRExA5gPXF47lpkJ\nXAYcMMnTDigfr3fpFOdLkqQe63vIAHYBtgVuazh+G7DbJM/Zrc3zJUlSjw3yjp8BZAfP3wHgoovW\nsnr15g/ceSfs0rT7o9pjnX69YX6vDRs2sLrxA1dX+Zn3np957/mZ99batWtrd3do53mDEDLuBB4G\ndm04PostRytqbm3zfIAnAnz4w29ov0JNy/z58/tdwtjxM+89P/Pe8zPviycCV7d6ct9DRmY+GBGr\ngIOAiwEiIsrvT53kad9r8vjB5fHJXAq8HrgJuH96VUuSNFZ2oAgYl7bzpCh6LPsrIg4HzgH+GlhJ\nsdrktcBTM/OOiDgX+GVmvr88/wDgSuB44N+AJeX9eZn50z78CJIkqUHfRzIAMvPCck+MkyimQa4B\nDs3MO8pT9gAeqjv/exGxBPhIefs58BcGDEmSBsdAjGRIkqTRMwhLWCVJ0ggyZEiSpK4Yi5DR7sXX\nND0R8fyIuDgifhURGyPiVf2uaZRFxPsiYmVE/C4ibouIiyLiT/td1yiLiCPLCzNuKG9XR8RL+l3X\nOCn/3W+MiFP6Xcuoiohl5Wdcf2ur93HkQ0a7F19TRzyKonn3aNrbUE3VPB84DXgW8GJgBvAfEfHI\nvlY12tYDf0txSYT5wBXAv0bE3L5WNSbKPxTfTvH/c3XXtRQLMnYrb89r58kj3/gZESuA72fmu8vv\ng+J/EKdm5if7WtwYiIiNwGGZeXG/axkXZYC+HViYmd/tdz3jIiJ+A/xNZp7V71pGWUQ8GlgFHAV8\nAPhhZh7X36pGU0Qso1i5Oa/qa4z0SEbFi69Jw+4xFCNId/W7kHEQEdtExGJgR6beEFCd8Tnga5l5\nRb8LGRNPKae+b4iI8yJiz3aePBD7ZHTRVBdf26f35UjdVY7UfRr4rvvGdFdE7EsRKnYA7gb+MjOv\n629Vo60Mc88A/rzftYyJFcCbgJ8BuwMnAFdFxL6ZeU8rLzDqIWMy7V58TRoWpwNPA57b70LGwHXA\nfhQjR68Bzo2IhQaN7oiIPSgC9MGZ+WC/6xkHmVm/hfi1EbES+AVwONDStOCoh4wqF1+ThlJEfBZ4\nGfD8zLyl3/WMusx8CLix/HZ1RCwA3k3RK6DOmw/8P8CqcsQOipHqhRHxLmD7HPUmwz7LzA0RcT3w\n5FafM9I9GWXarV18Ddjs4mstX0VOGnRlwPgL4IWZeXO/6xlT2wDb97uIEXYZ8GcU0yX7lbf/As4D\n9jNgdF/ZdPskoOU/YkZ9JAPgFOCc8kqvtYuv7Qic3c+iRllEPIoi6db+2pgTEfsBd2Xm+v5VNpoi\n4nSKiwS+CrgnImojdxsy0ysOd0FEfAT4OsVKtZ0orvB8IHBIP+saZWUPwGZ9RhFxD/CbzFzbn6pG\nW0R8CvgaxRTJE4ATKa4jNtHqa4x8yGjh4mvqvD8HvkXR95IU+5RAcaXdt/SrqBF2JMXn/O2G428G\nzu15NeNhV4rPdndgA/Aj4BBXPPScoxfdtQfwRWBn4A7gu8CzM/M3rb7AyO+TIUmS+mOkezIkSVL/\nGDIkSVJXGDIkSVJXGDIkSVJXGDIkSVJXGDIkSVJXGDIkSVJXGDIkSVJXGDIkSVJXGDIkSVJXGDIk\nSVJXGDIk9UVE7BIRt0TE8XXHDoiIP0TEC/tZm6TO8AJpkvomIl4KfBU4APgZsAa4KDPf09fCJHWE\nIUNSX0XEacDBwH8B+wL7Z+aD/a1KUicYMiT1VUTsAFwL7AHMy8yf9rkkSR1iT4akfnsS8HiK/x/N\n7nMtkjrIkQxJfRMRM4CVwA8pejKOA/bNzDv6WpikjjBkSOqbiPgU8Grg6cC9wLeB32XmK/tZl6TO\ncLpEUl9ExIHAscAbMvOeLP7ieSPwvIj46/5WJ6kTHMmQJEld4UiGJEnqCkOGJEnqCkOGJEnqCkOG\nJEnqCkOGJEnqCkOGJEnqCkOGJEnqCkOGJEnqCkOGJEnqCkOGJEnqCkOGJEnqCkOGJEnqiv8fHH2P\ndQopiUwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb48bd60cf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def pareto_dist(x, xm, alpha):\n",
    "    '''\n",
    "    x: iterable of x values\n",
    "    xm: parameter x_m\n",
    "    alpha: parameter alpha\n",
    "    '''\n",
    "    pdf = []\n",
    "    for point in x:\n",
    "        if point < xm:\n",
    "            pdf.append(0)\n",
    "        else:\n",
    "            pdf.append(alpha*(xm**alpha)/point**(alpha+1))\n",
    "    return pdf\n",
    "\n",
    "# lambda value\n",
    "xm = 1\n",
    "alpha = 1\n",
    "_ = plt.plot(x, pareto_dist(x, xm, alpha), '.', markersize = 1)\n",
    "_ = plt.ylabel(\"Probabilty mass at x\")\n",
    "_ = plt.xlabel(\"x\")\n",
    "_ = plt.title(\"Pareto Distribution\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAFyCAYAAACKkcLXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmcFNW5//HPM4jsIoosV0EGVKLGoKDEGBm2uFzDNYYY\nAX+5irhEk4Bwo4maRBjRRJMrJpioxETFaABFkxiXuAtqVHKHuCW4xBkUFUUU0ABucH5/nKqZ6pru\nme6e7unt+369+gVTXV11urq66ulznnOOOecQERGRylZV6AKIiIhI4SkgEBEREQUEIiIiooBARERE\nUEAgIiIiKCAQERERFBCIiIgICghEREQEBQQiIiKCAoKKYGbbzezCQpcjysweMbNnC10OyY6ZzTGz\n7YUuR7kysz2D7+1JWbx2avDagVnu+wYza4gta5driJmNDvZVE1nWbteKthz3clARAYGZDTazBWb2\nipltNbNNZvaYmc0ws86FLl87cMGjXZlZfzObbWafS1GmimVm55vZVwpdjjYoyDlVbsxsipmdneLp\nbI9vWz+bZK/PeJutvLfW9t/S322Wp+Ne8nYodAHyzcyOAW4FPgRuBJ4HdgQOB34K7AecWbACto8u\nwKcF2O9/ALOBBkC1AYkuwJ+Xfyp0QaSgTgT2B34RXeice9XMugCfFKRUzWVzDUn63lrinFtmZl2c\ncx9nuK9Mlcpxb1dlHRCY2SBgMf6GNM45ty7y9NVm9iPgywUoWt6ZmQE7Ouc+aocvV8piFGi/7crM\nOgBVzrmKvIhIfhTwe9tMvstiZp2Aj51X0Pdd6P0XlHOubB/A1cA24PNprt8B+BHwL3yNQgNwMf7G\nGl1vNXAHMBr4G7AF/wt4dPD8xODvrcD/AQfGXn8D8AFQDdwL/Bt4A/hRkjKdAzwOrA/283/A15Ks\ntx2Yj498nwc+Ao6NPHdhZN05wbIhQVk2ABuB64DOse12Drb7DvA+8Ef8L/+EbSYpz+hgnW3Bv+H/\nTwqefzg4RvsG/98MvA6cm2RbOwK1wMvB5/IacFn8c0lRjkeC/QwPjuMWoB74Zmy9jsBFwfHdGHwm\ny4ExsfX2DN7L/wBnB+fKJ8Dn0i1rkuOyHbgu8vxBwD3ApuA8eYD0z+FMz5evAM8FZX0eOCrJuofj\nz/Otwfs6A1/zsy3NMn0e+EtwXDcHn8lhsXNsVfDoFFneC1gLPJbld6crcHnwGXwIvAB8t43H4j/w\n35O3IutNS3Hufx34AbAmOHYPAEMi6z2c5Fyoj51nJ0XWPwC4Hngl2N5a4LfALrH9nxxsc2Aan81x\nwXvYiv+eHBfsoz7JMYpeQ7oDP8dfIz8E3gbuI7jWtfLexgR/T8JfX9fgax92ihy7mthxSuc7PDV4\n7cAUn0dNNsc9WD4OeDQ43zbgr4Ofia0zh/Svq0cE29uAP59fAC5J5/uUz0dZ1xAAE/Af9FNprv9b\n4CTgFuB/8ReyC/A3ra9F1nPA3sDNwALgd8C5wB1mdhZwCfAr/C/kC4AlwNDY66vwF8kngtceDdSa\nWQfn3JzIujPw1co34W82k4FbzGyCc+6eWPnH4y9Cv8LfEFaneJ9hG9kt+C/Wefgv22n4L/b5kXUX\nAsfjm1uewn+57qL1drZVwIX4m+wC/MkP8NfIOrvgb3y342tyjgcuNbNnnXP3QmNNx5+Bw4LtvIC/\nMM7CfwYTWymHC/ZzV/B+fw+cgK8h+sg5d0Ow3k7ANGAR8GugB3Aq8BczG+mcizd5TAM6BWX6CHgv\ng7J+A3+uPRXsC/xFHjPbHx+IbAIuxV8ovwk8YmY1zrm/tfJ+MzlfRgVlugp/UZoBLDWzPZ1z7wXl\n+Sz+xrsO/3l2xF/41pEGMxsH3I0PTObgL5inAA+Z2eHOuf9zzn1oZifjL/aX4IMagnL1wN/gQpl8\nd/6MP19/CzwNHAX8zMz+wzn33SyORR/8Z7YNH0CsB/4T+I2ZdXfOzY9t87xg3Z8BPYHv4z+XLwTP\nXxws3x2Yib9e/LuFw3kEPhAKA5L98efGfpFtps3MjgSW4gOC84Bd8cHA62m8fAH+eF2J/67vCnwR\nf618mpbfW3jt+BH+u/O/+O/Sx7Hno9L5DreU5xBdntFxN7Mv4c/hV/CBcBf8+fGYmQ13zr0W20eL\n11Uz2w9/bj4dOQZ74a8bhVXoiCRfD/yFZDtwe5rrfy5Y/5rY8p/iv9SjI8saiNU84L+s2/En1h6R\n5acH60Yj3uuDZVfE9vVnfKS+S2RZp9g6HfDR8v2x5dvxv1SHJnlv8eh+drDs17H1bgPWRf4+KFjv\nf2PrXReUP2UNQbDeCJJE28FzDwfbODGyrCP+V88tkWXfCN7XF2KvPyN4/aGtlCHcz9mx/awM9tUh\nWGbADrHX7hSsc21kWfgLYgPNf5mlXVb8Tee6JOX9Q3AO7BlZ1g8fIDycxnmcyfmyFRgUWXZAsPxb\nsfJsBnaPLBsavM9WawiAF4G74mXEX1z/Elt+SbDdL+KDw+3Ad2LrpPXdwf/a3w6cF1tvCT7Iqs7i\nWPwGf7PcObbN3wPvhceepl+kz4fnV7B8elD2/WLlTvg1HjvPojUEnZKsNynY5hcjy9KqIQD+Hryf\n7pFl44n8Yo4do+g1ZAMwv5Xtp3pv4fF5mea1r6Npfr1M9zuc9H2n2GYmx/3vwX56xs6PT4HrI8vS\nva6eHZSnV2vfn/Z+lHMvg52Cfz9Ic/1j8BHeFbHll+NvFvFcg3+6xJqH8P8POudejy03YHCSff4q\n9vcv8b/qvhQucM59FP7fzHbGV6M+io884x5xzr2YZHkyDh/lRz0K7Gpm3YO/jw7Wuzq23pXkJj9g\ns3Pu940F8m3wT5F4rI7H/wJ5ycx2DR/4i4QBY9PYz6c0/RIP97MA6IMPWnDep+BrJcysF/6z+D+S\nH+ulLvjlmKuymlkVPrD8g3Pu1Uh538LfdEZFPpukMjxf7nfOrY689jl8s9DgWHn+6Jx7I7Lei/ha\ngxaZ2YH4mpFFsePRA3gQqIm9ZA7+Jnoj/rvxsHPulyk2n+y704mm784x+M/9yth68/A1DP8ZW97i\nsQhMxN9IOsTez334X5zxY3ydc25b5O9HSX0taFXss+0U7Du8viT7fFMys37AMOAG51zjr2Pn3IPA\nP9PYxEZgpJn1z2S/MTe49NvrW/0O50PkOF3vnNsU2f9zwP348ywqnevqxuDfrwa1ikWjnAOC94N/\ne6S5fhgZ/iu60Dn3Nv4D3DO2/mux9cL9xavbwpOoV2z5dny1UtRL+C93477MbIKZPWFmW/G/QtYB\nZ+EvQHGrkyxryWuxvzfEyhoek4bYev8iN9YkWbaBxGO1N75q9J3Y40X8l69PGvt50zm3NbYs2bE+\n2cyewbeJvos/1l8m/WPd1rLuhm/3finJc6uC8g5oaQMZni+tHf+wPC8nWS+dwHPv4N8bSTwe6/DV\nqDuaWWO5gov8afhq8e74ZplkUn13oOnzHIj/3DfH1lsVWy/U4rEws92AnfG1PfHP97pg/fjnG99m\n/PuVETPrZWa/MLO38DUa7+CPgyP559uS8P0n+y6n89l+D/8reY2ZPRV0L67OsAyrM1g3re9wHoTb\nTvWd7B30Sohq7bq6BN88di3wtpktMrOvF0NwULY5BM65D8zsTfxJm47ww2itbTy0LcPl6XzYCeuY\n2Sh8e/Aj+Iv6WnyV6jRgSpLXx78wrcm2rOkeo1zsvwqf6DWL5OVKdiFPR/xYfwNfHX07vploXVC+\nC0j+iy7ZsW5rWdt0QcjifGnt+Lf0nUinrOEPju8Cz6RYJ952e3Twb2d8QPEq6YmXJ1X5Up27rR2L\n8L3chM+rSSaeZ9KWa0EytwKH4s/PZ/DHrgpfW5Ppj7s2fbbOuVvNbDnwVeBIfN7H983sqy7I/0lD\npteruHg5U322HXK4j3S0+Lk75z4EasxsLP4Hx9H4pp8HzexIF7QrFELZBgSBO4HTzezzrvXEwtX4\nL9XeRCLkIJFoZ9K/MKWrCn+jiUbo+0TKAj6RcSs+27mxD7CZnZrjsqTyKr6c1QRJb4F9kq/eTC5O\n7FfwGfwPt2Eb/xH0bY5egPbBly/8XL8GvOKcOz76QjO7KE9lTXZs1uEzqIcmeW7f4DUtBRUTye35\nsi7YXrLPO1kZ48Jz5gPn3EOtrWxmB+CTrK4DDsQn6x3gnIs3+6Xz3VkNjDOzbrFagv2CfzP9Pr+D\nb37skM57yUBa35Gg+WccvjfFJZHle2W539XBv8k+27S+30Ht6TXANWbWG9/W/gOampNyeWNL5zsc\n/hLfmcRf6YOSbC/dsq0O/k12vn8GWJ+k5iItwXXiYeAcMzsfn+w4Fsjl+ZWRcm4yAB9Jb8FfWJpV\n15rZEDObEfx5Nz6Cmxlb7bv4k+euPJTvO0n+/pimE+LTYN+NgVswtkJ7jXB3L/6YfCu2fDrpfaHC\nC/HObSjDLcAeZnZ6/Akz62xmXdPYxg5EBp8ys4747Ox3gLpg8TZi78nMPk9m2duZlHUzsePinNuO\nb4/+ikWGnTWzvvhf+Muj7b1JhO8hJ+dLUJ57gePMbI/INvfF/ypsTR0+KDjHzLrFnwxuIuH/d8D/\n8n4dn3Q1FZ9MGc/pCbX23bkbfxzi683CNznEe1y0KDgWtwFfC3qCJIi+lwxtJr3q/vBXZ/yaPYss\nbrxBXsrTwMlm1tisamZH0BQ0JWVmVWa2U3SZc2498CY+jyOU7ntLRzrf4Vfw16vosMdV+GaeuLTK\nFjtOje856H1zJFncF4L8pLhn8GXvlOS5dlPWNQTOuXozOxHfpW2VmUVHKjwM30Xv+mDdZ81sIXBG\n8IEtw3c7PAnfU2FZjov3EXB0sM8n8ckp/4nvi/pusM6d+P7u95rZ74G++Jvzy/heEXnlnFtpZrcB\nM4ML3pP4jN2wbbi1C9Er+PyLM83s3/gv4ZPRhLk0/I6mLkZj8W1vHfC/mL+O/1KubGUbbwLfC9o4\nX8R3xfsccHok6etOYKKZ/RH/JR+Mv+D8A9+eneuy1gFfMrNZQfkanHMrgB/iE+MeN7Or8DeCM/Dn\n7Pda2X8+zpfZ+CrNx4LydMTfZP9BK81xzjlnZqfhb87/MLPr8WMG7I7/JbSJpmDlR0EZxwW/6J8P\namcuNrPbXGKXyVa/O865O8zsIeASMxtMU7fD/8L3UIjnxaTjPHwf+qfM7Fp88t0u+KS2cUA2QUEd\ncIKZXY4f6+Hfzrk74ysFTaDL8efxjvjjeCS+9i7bJojz8efM42Z2Hb7r4Hfw18iWzvkewOtmtpSm\nposjgIPx519G7y2JZO9nLa18h51z/zSzJ/Bdl3fF59BMJvkP30zKdi7+HH7SzH6Lz6v5Dr5GojaN\n9xN3ofm5Gu7C1270xTfxvQY8lsX2cqfQ3Rza44EfKOIamgb02Ijv630W0DGyXhX+ghwOTLQamBtd\nJ1ivHvhTkv1sA34RW7ZnsHxWZNn1+KTHQfj+1B/gbwrJBleZiu/PvgV/ET6JJIPCJNt37LkfRf6e\nHSxrdUATEgcm2oTvt7wX/ldWs0GEkux7Ar5d/SOaD0z0TJL1r8dX3UeXdcC3UT4bHIf1wAp89WT3\nVvb/cPC6g/A36M3B53dmknW/HzwXDujzn/HyJPs8sykrvrrzYfzFdBuJAxMNw1+AwoGJ7gdGpnmu\nt+l8Cd7/b2PLDg/eQzgw0enJttlCmT6Hb/8OmyDq8eM9jAmePyg4P+JdCavwWfRrgJ2y+O50xfdx\nX0PTwETNPrcMj0Vv/PdhdbDNN/C1OtMi64Td3CamuBacFCvj7/BJrNtIHCAnvm5//PfvXfzNbhH+\nZhL/fiftfpfiswkHJtqC/55+heTfwcZ94IPCS/HB7cbg81gJnJHk+Cd7b0mPT+y5eLfDZ0jvOzwI\nX6u1JTgvLsIHa/Ftpn3cg+Vj8feMcGCiPxDr4k2a19VgW7fTNGDVmqAsQ+Lvp70fFhRQ2lHwS+lr\nzrmdWl25CAXdyVYC/885t6jQ5WmJmT0M7Oqcy3uNiuRfqX93RIpZxjkEZjbKzO4wszfMTxN5bBqv\nGWNmdWb2oZm9ZH5UMikB5scYj5uJj3iXt3NxREQkT7LJIeiGb4+7Dp9k06IgqelO/JCgJ+LbR39j\nZm865+7PYv/Svr5nZiPwXdk+xbfXHgUscJHBakREpLRlHBA45/6Cb7sLx5lvzVn49pkwIepFMzsc\nnx1byQFBqbTVPIFPGPohPtHoNXxb2Y8LWagMlcqxlvTo8xTJgzblEJjZduA459wdLayzDKhzzv1P\nZNlUfAJRViN2iYiISG61R7fDfviZnqLeBnYys04uMj53KOgychRNmbwiIiKSns4EPS5cUzf2VhVq\nHILWhgk+Cj+1sIiIiGTn/+EnRktLewQEb+H7ykb1Ad53qWe6Wg1w0003se++++axaBJ64/03OPbk\nY7lj4R3svtPuhS5OxZg1axZXXJFqMD7JBx3z9qdj3r5WrVrFN77xDchwwrv2CAieoPlUo0cGy1P5\nEGCnnfZl+PCMZvWULA1nOOP3G89/jfmvQhelovTs2VPneDvTMW9/OuYFk1GTezbjEHQzs2HB4DQA\ng4O/BwTP/yQYUjR0DTDEzC4zs6Fm9i38vPHzWtvX7vqhKiIi0i6ymdzoYPysVnX4HIDL8aPWhWM6\n9yMyZ7tzbjV+iscv4ccvmAWc6px7oLUdrWxthHrJmYYNDTxY/yANG7IZ4l1EREpdNuMQLKOFQMI5\nd0qK14zIdF+nnw7jx0N1daavlExV96pm/ODxVPfSwRYRqURFPf3xtdcqGGhPp550qmoI2tmUKVMK\nXYSKo2Pe/nTMS0NRTm5kZsOBOqhj2bLh1NS0+hLJgYYNDQyeP5j6GfWqKRARKVErV65kxIgRACOc\nc2k3vhdqHIK0XHstjB4N9fWqKWgP1b2qFQyIZOC1115j/fr1hS6GVJjevXszcODAnG+3qAOC4cNh\n2TIFAyJSfF577TX23XdftmzZUuiiSIXp2rUrq1atynlQUNQBwRtvwLHH+qBAzQb5pyYDkfStX7+e\nLVu2aAA1aVfhoEPr16+vrIBg9919MKBmg/ahJgORzO27rwZQk/JQ1L0MwNcMqNmg/VT3qlZPAxGR\nClT0AUFDg68hWL680CWpDGGzgYICEZHKUvQBQXV1U7NBg+5ReadmAxGRylT0AQGo2UBERCTfSiIg\nULNB+1GTgYhI9saMGcPYsWOzem1VVRUXXXRRjkuUwf4LtucMqNmg/ajJQEQAFi5cSFVVVdJHhw4d\nWLFiRaGLWDCrVq2itraW1157rdlzZkZVVUncWpsp6m6HUWo2EBFpX2bG3LlzGTRoULPn9tprr/Yv\nUJH45z//SW1tLWPHjm02FsD9999foFK1XckEBGGzgQYpyi8NTiQiUUcffbTGWYhxzmFmSZ/bYYeS\nua02UzL1Gmo2aB9qMhCRdM2ePZsOHTrw8MMPJyw//fTT6dSpE8899xwAy5Yto6qqiltuuYULLriA\n/v370717d77yla/w+uuvN9vurbfeysEHH0zXrl3Zbbfd+O///m/efPPNhHWmTp1Kjx49ePPNNznu\nuOPo0aMHffr04dxzzyU+aZ9zjp///Od89rOfpUuXLvTr148zzzyTjRs3Jqw3aNAgjj32WB5//HE+\n//nP06VLF4YMGcLvfve7xnUWLlzICSecAPh8gbAJZXmQ5DZmzBjGjRvXuP4nn3zChRdeyMEHH8zO\nO+9M9+7dqamp4ZFHHsnwaOdfyQQEAAMGFLoElUGDE4lIaNOmTbz77rsJj/feew+AH/3oRxx44IGc\neuqpbN68GYB7772X3/72t8yZM4cDDjggYVuXXHIJ99xzD+eddx5nn302999/P0cccQQfffRR4zo3\n3HADkyZNomPHjlx66aWcccYZ3H777YwaNYr333+/cT0zY/v27Rx11FHstttuXH755YwZM4Z58+bx\n61//OmG/Z5xxBt///vcZNWoU8+fPZ9q0adx8880cffTRbNu2LWGbL7/8Ml//+tc58sgjmTdvHrvs\nsgunnHIKq1atAqCmpoYZM2YA8MMf/pCbbrqJ3/3ud43DV8drDt5//32uu+46xo4dy09/+lNqa2tZ\nv349Rx99NM8++2ybPpucc84V3QMYDri6ujoXV1/vH5I/9e/VO+bg6t/TgRZJpa6uzqW6TpWDG264\nwZlZ0keXLl0a13v++eddp06d3BlnnOE2btzodt99d/f5z3/ebdu2rXGdRx55xJmZGzBggNu8eXPj\n8ltvvdWZmbvyyiudc8598sknrm/fvm7YsGHuo48+alzvrrvucmbm5syZ07hs6tSprqqqyl1yySUJ\n5R4+fLg75JBDGv9+9NFHnZm5xYsXJ6x33333OTNzixYtalw2aNAgV1VV5R5//PHGZe+8847r3Lmz\nO/fccxuXLV261FVVVblly5Y1O25jxoxxY8eObfx7+/bt7pNPPklYZ9OmTa5fv37utNNOS1huZq62\ntrbZNqPSOe/CdYDhLoN7b0nVEIQGD1azQT5V96pm2cnL1GwgkmP5vG7lY9tmxtVXX80DDzyQ8Ljn\nnnsa19l///2pra3l2muv5aijjuK9995r7KEQd/LJJ9O1a9fGv48//nj69+/P3XffDcDf/vY31q1b\nx7e+9S123HHHxvWOOeYYPvOZz3DXXXc12+Y3v/nNhL9HjRpFfX19499Lly5l5513Zvz48Qm1HAcd\ndBDdu3dv1tyx3377cdhhhzX+3bt3b4YOHZqwzUyYWWNegXOODRs28PHHH3PwwQezcuXKrLaZL6Wb\n/SB507ChgdELRyuXQCSHGhr8j5l8TNSWz20fcsghrSYVnnvuuSxevJi//e1v/PjHP2bo0KFJ10vW\nM2Gvvfbi1VdfBfyU0mbGPvvs02y9z3zmMzz++OMJyzp37syuu+6asKxXr15s2LCh8e+XX36ZjRs3\n0qdPn2bbNDPWrVuXsCzZDILxbWZq4cKFzJs3jxdeeIFPPvmkcfngwYOz3mY+lFxAECYXqvth/iix\nUCT3qqvzN2trPredjldeeYWXX34ZoDGRMF0ukgAY/X86OnTo0Oo627dvp2/fvvz+979Puv3ddtst\nrW1mWrbQTTfdxCmnnMLEiRP53ve+R58+fejQoQM//vGPs651yJeSCwjU/bB9hImFCgpEciefN+xC\nBQPOOaZOnUrPnj2ZNWsWl1xyCccffzzHHXdcs3XDoCHqlVdeYdiwYYDP8nfO8eKLLzJmzJiE9V58\n8UX23HPPjMs3ZMgQHnzwQQ477DA6deqU8euTSdXlMJnbbruNIUOGsHTp0oTlF154YU7Kkksll0Og\n7oftQ0MYi0g6Lr/8cp588kmuvfZaLrroIr74xS9y1llnNfZEiLrxxhv597//3fj3rbfeytq1aznm\nmGMAOPjgg+nTpw/XXHNNQtX6Pffcw6pVq5gwYULG5TvhhBP49NNPkw4JvG3bNjZt2pTxNrt164Zz\nrlm3xWQ6dOjQLIB46qmneOKJJzLeb76VXA0BaNTC9qBmAxFxznH33Xc3drmLOuyww/jwww+58MIL\nOeWUUxpv6tdffz0HHnggZ511FkuWLEl4zS677MLhhx/OKaecwltvvcUvfvEL9tlnH0477TTAD+pz\n2WWXMW3aNGpqapgyZQpvvfUW8+fPZ/DgwcycOTPj91BTU8M3v/lNLr30Up5++mmOPPJIOnbsyEsv\nvcTSpUuZP38+EydOzGibBx54IB06dOCyyy5j48aNdOrUifHjx9O7d+9m606YMIHbb7+d4447ji9/\n+cvU19ezYMEC9t9//4TgqBiUZECgZgMRkfwzM2bPnp30ud/85jdcc8019OnThyuuuKJx+V577cVP\nfvITZs6cydKlSzn++OMbt3XBBRfw7LPPcumll/LBBx9wxBFH8Ktf/YrOnTs3vv7kk0+mW7duXHrp\npZx33nl069aNr33ta1x66aXstNNOzcqXqtxRV199NQcffDALFizgBz/4ATvssAODBg3ipJNO4otf\n/GLC69LZZt++fVmwYAE/+clPOO2009i2bRsPP/wwNcENKbru1KlTefvtt1mwYAH33Xcf++23Hzff\nfDO33HJL42BG6ey/PVi2iRL5ZGbDgbq6urqU2a3Ll/ugoJCJNOVMQxiLtGzlypWMGDGClq5T4i1b\ntoyxY8eydOnSjH+NS6J0zrtwHWCEcy7tvo0ll0MQ0qiF+RU2GYiISGUo2YAg7GYj+aXEQhGRylCy\nAUFIoxbmjxILRSRXCtk2LukpyaRCEREpHaNHj06YREiKU0nXEKjZIL80FoGISOUo6YAgpGaD/FCT\ngYhI5SiLgABgzZpCl6A8hUMYi4hIeSv5gEBDGeeXmg1ERCpDyQcEoKGM86m6VzXLTl6mZgMRkTJX\nFgFBOJRxbBRIyYGGDQ2MXjhaNQQiImWuLAICNRvkjxILRUQqQ1kEBKChjPNJiYUiIuWvbAICjUmQ\nP0osFKk8CxcupKqqqvHRpUsXhg4dyvTp01m3bl27l2fr1q3U1tY2myFQcqfsRiocPFgzIOaaEgtF\nKpOZMXfuXAYNGsSHH37IY489xtVXX80999zD888/nzBtcb5t2bKF2tpazKxxmmHJrbKpIZD8UWKh\nSOU6+uijOfHEE5k2bRrXXXcdM2fOpKGhgT/96U9t2u5HH32Ecy7t9TNZV7KjgEBapcRCEQmNGzcO\n5xwNDQ1s2LCBc845h8997nP06NGDnj17cswxx/Dss88mvGbZsmVUVVWxZMkSfvjDHzJgwAC6devG\nBx98AMCmTZuYOXMmAwcOpHPnzuy999789Kc/bQwCXn31Vfr06YOZMWfOnMZmjIsuuqhxHw899BCj\nRo2ie/fu9OrVi+OOO44XXnih/Q5MGSirJgPlEeRPmFiooECksv3rX/8CYNddd6W+vp477riDr3/9\n61RXV/P222+zYMECxowZwz//+U/69euX8Nq5c+fSqVMnzjnnHD766CN23HFHtm7dSk1NDW+++SZn\nnXUWAwYM4K9//Svnn38+b731FvPmzWO33Xbjmmuu4cwzz2TixIlMnDgRgM997nMAPPDAAxxzzDEM\nGTKE2tp2ZKfUAAAgAElEQVRatm7dyvz58zn88MNZuXIlAwcObN+DVKqcc0X3AIYDrq6uzmWqvt45\n8P9K7tS/V++Yg6t/TwdWxDnn6urqXLbXqVJwww03uKqqKvfQQw+59evXu9dff90tXrzY9e7d23Xr\n1s29+eab7uOPP272uldffdV17tzZXXzxxY3LHnnkEWdmbq+99nIfffRRwvpz5851PXr0cK+88krC\n8vPPP9917NjRvf76684559avX+/MzNXW1jbb54EHHuj69evnNm7c2Ljs2WefdR06dHBTp05t03Eo\nNumcd+E6wHCXwb23bJsMNLdBbimxUKTt8pmHk49tO+cYP348u+22GwMGDODEE09kp5124o9//CP9\n+/enY8eOjetu376d9957j65duzJ06FBWrlzZbHtTp05lxx13TFi2dOlSRo0aRc+ePXn33XcbH+PH\nj+fTTz9ttVfBW2+9xTPPPMMpp5xCz549G5cfcMABHHHEEdx9991tPAqVo6yaDCBxkCL1NsidMLFQ\nuQQi2Qm77+bjO5SvbZsZV111FXvvvTc77LADffv2ZejQoY3PO+f4+c9/ztVXX01DQwPbtm1rfF3v\n3r2bbW/QoEHNlr388ss899xz7Lbbbkn331oXx1dffRWAffbZp9lz++67L/fddx9bt26lS5cuLW5H\nyjAgAM1tkA+qIRBpm3wm5+Zz24cccgjDhw9P+twll1zChRdeyKmnnsrFF1/MLrvsQlVVFWeffTbb\nt29vtn6ym/L27ds54ogj+P73v5+0J0GyG31UstdIdsoyIAjnNli2zAcH0naqIRBpu3x+dwrxvbzt\nttsYN24c1157bcLyjRs3Jv3Fn8yQIUP497//zdixY1tcz8ySLg9rHV588cVmz73wwgv07t1btQNp\nKsscAs1tkHvqeigicR06dGj2C/3WW2/ljTfeSHsbJ5xwAk888QT33Xdfs+c2bdrU2AzRtWtXwAcb\nUf369ePAAw9k4cKFvP/++43Ln3/+ee677z6+/OUvp12WSleWNQTQNLfBmjVqOsgVdT0UqSytVcdP\nmDCBuXPnMm3aNA477DCee+45br75ZoYMGZL2Ps4991zuuOMOJkyYwNSpUxkxYgSbN2/m2Wef5fbb\nb2f16tXssssudO7cmf32248lS5awzz770KtXLz772c+y//7787Of/YxjjjmGQw89lFNPPZUtW7bw\ny1/+kl69ejF79uy2HoaKUZY1BKBagnzQnAYilSVVNX3oggsu4Lvf/S733XcfM2fO5Omnn+buu+9m\nwIABzV6baltdunRh+fLlfO9732PZsmXMnDmTyy67jFdeeYWLLroooefAb3/7W3bffXdmzZrFiSee\nyG233QbA+PHj+ctf/kLv3r2ZPXs28+bN47DDDuOxxx5jzz33bONRqBxWjAkZZjYcqKurq0uZzJKu\n5cuVR5BLy1cvp2aQDqjIypUrGTFiBLm4TomkK53zLlwHGOGca97/M4WsagjM7Ntm1mBmW83sSTM7\npJX1Z5rZC2a2xcxeM7N5ZtYpm31nIkwu1ORYuaE5DUREylfGAYGZTQIuB2YDBwHPAPeaWfNOp379\nE4GfBOt/BpgGTAIuybLMaVOzQW6FiYUiIlJ+sqkhmAUscM7d6Jx7ATgT2IK/0SfzBeAx59wS59xr\nzrkHgEXAyKxKnKFocqHkhvIIRETKT0YBgZl1BEYAD4bLnE9CeAB/40/mr8CIsFnBzAYDxwB3ZVPg\nTKmWILfU/VBEpDxl2u2wN9ABeDu2/G1gaPPVwTm3KGhOeMx8mmkH4Brn3GWZFjZbGrlQRESkZbka\nh8DwMys1f8JsDHABvmlhBbAXMN/M1jrnLm5po7NmzUrocgIwZcoUpkyZklHhNHJh7uRzPHYREcnM\nokWLWLRoUcKyTZs2ZbWtTAOC9cA2oG9seR+a1xqELgJudM5dH/z9DzPrDiwAWgwIrrjiipx059GE\nR7mjxEIRkeKR7EdypNthRjLKIXDOfQLUAePDZUEzwHh8rkAyXYH4LBfbg5e2POpFDoXJhZIbSiwU\nESkv2TQZzAMWmlkdvglgFv6mfwOAmd0IvO6cuyBY/8/ALDN7GngK2Btfa/AnV4yjIkmrNPOhSJNV\nq1YVughSQfJ5vmUcEDjnbgmSBC/CNx08DRzlnHsnWGUP4NPIS+biawTmArsD7wB3AD9sQ7kzVl3t\nmwuk7TTzoQj07t2brl278o1vfKPQRZEK07VrV3r3Tjr0T5tklVTonLsKuCrFc+Nif4fBwNxs9pVr\ngwcrj6CtlEcgAgMHDmTVqlWsX7++0EWRCtO7d28GDhyY8+2W7WyHLdEMiLmh3gZS6QYOHJiXC7NI\nIZTtbIfJaJCi3FEegYhIeamogAA0lHGuaKIjEZHyUnEBgWoJckN5BCIi5aXiAgLQUMa5pPEIRETK\nQ0UGBOFQxsuXF7okpU15BCIi5aMiAwI1G+SG8ghERMpHRQYEoOTCXNBUyCIi5aNiAwLVEuRGda9q\n1RCIiJSBig0IQBMe5UI4HbKCAhGR0lbRAYHmN2g7JRaKiJSHig4IQoMHq9kgW0osFBEpDwoIpE00\nQJGISHlQQBBQb4O2UR6BiEhpq/iAQL0NREREFBAAGpOgrZRYKCJS+hQQoFqCtlJioYhI6VNAENCY\nBNlTYqGISOlTQBDQmARtp8RCEZHSpYAgRmMSZEd5BCIipU0BQRJKLsyc8ghEREqbAoIIJRdmT3kE\nIiKlTQFBjLogto3yCERESpMCghjVEmRPeQQiIqVLAUESqiXIjvIIRERKlwKCJFRLkB3lEYiIlC4F\nBCmoliB7yiMQESk9CghSUC1BdpRHICJSmhQQtEDDGWcuzCNYvnp5oYsiIiIZUEAgORXWECi5UESk\ntCggaIHmN8hOzaAaNRuIiJQYBQRp0PwGmVH3QxGR0qOAIE3qbZA+dT8UESk9Cghaod4G2VP3QxGR\n0qGAIA0akyBz6n4oIlJaFBCkQbUEmVP3QxGR0qKAIE2qJciMuh+KiJQWBQRpUi1B5tT9UESkdCgg\nyIBqCTKj7ociIqVDAUEGVEuQGXU/FBEpHQoIMqRagsyp+6GISPFTQJAh1RKIiEg5UkCQBdUSpE/j\nEYiIlAYFBFlQLUH6NB6BiEhpUECQpbCWQFqm8QhEREqDAoI2UrNB6zQegYhI8VNAkCU1G6RP4xGI\niBQ/BQRtoOTC9Gg8AhGR4qeAoA1US5AZjUcgIlK8FBC0kWoJRESkHCggaCPVEqRHzQYiIsUtq4DA\nzL5tZg1mttXMnjSzQ1pZv6eZ/crM3gxe84KZHZ1dkYuPagnSp2YDEZHilHFAYGaTgMuB2cBBwDPA\nvWbWO8X6HYEHgIHARGAocDrwRpZlLjqqJRARkVKXTQ3BLGCBc+5G59wLwJnAFmBaivVPBXYGjnPO\nPemce80596hz7rnsilycVEvQOjUbiIgUr4wCguDX/gjgwXCZc87hawC+kOJl/wU8AVxlZm+Z2XNm\ndr6ZlVX+QrSWYLlG6W2Rmg1ERIpPpjfl3kAH4O3Y8reBfileMxj4erCv/wTmAt8FLshw30UvrCVQ\n04GIiJSaHXK0HQNciueq8AHDGUFtwt/NbHfgHODiljY6a9YsevbsmbBsypQpTJkype0lzoNoLYEk\np2YDEZHcWbRoEYsWLUpYtmnTpqy2lWlAsB7YBvSNLe9D81qD0Frg4yAYCK0C+pnZDs65T1Pt7Ior\nrmD48OEZFrGworkE1Rq6P6XB8wdTP6Ne8xuIiLRBsh/JK1euZMSIERlvK6MmA+fcJ0AdMD5cZmYW\n/P3XFC97HNgrtmwosLalYKBUqcdB+tZsUgamiEixyCaxbx5whpmdZGafAa4BugI3AJjZjWb248j6\nVwO7mtkvzGxvM/sycD7wy7YVvXipx0HLNCWyiEjxyTiHwDl3SzDmwEX4poOngaOcc+8Eq+wBfBpZ\n/3UzOxK4Aj9mwRvB/3/axrIXrWgtQX29mg6S0ZTIIiLFJaukQufcVcBVKZ4bl2TZU8Bh2eyrVCmX\noGXhlMjLTl5GzaCaQhdHRKTildVYAMVEuQQtU7OBiEhxUUCQR8olaJmaDUREiocCgjzS6IUtC5sN\nVEMgIlJ4CgjyTKMXtk7dD0VECk8BQZ6FtQSgpoM45RGIiBQPBQTtoKZGCYapKI9ARKQ4KCBoJ0ow\nTC7MI1i+WkkWIiKFpICgnagbYnJqNhARKQ4KCNpRWEsgiQb09AdGyYUiIoWjgKAAVqwodAmKi2oJ\nREQKTwFBO6quhsWLYfJkWLKk0KUpLkouFBEpLAUE7WzSpKagQLkETZRcKCJSWAoICmDkSP+vehw0\nUbOBiEhhKSAoAPU4SE7JhSIihaOAoEDCHgdKMGyiWgIRkcJRQFAgSjBMTsmFIiKFoYCggJRg2JyS\nC0VECkMBQYEpwTCRmg1ERApDAUGBKcGwuTC5UERE2o8CgiKgBMPk1NtARKT9KCAoAtEEw+VqOlez\ngYhIASggKBJhLoF4ajYQEWlfCgiKhHIJElX3qqZ+Rn2hiyEiUjEUEBQR5RI0N3j+YDUbiIi0AwUE\nRUSDFSWn5EIRkfxTQFBkNFhREyUXioi0HwUERShMMFTTgSY8EhFpLwoIipCaDpqolkBEpH0oIChS\najpooloCEZH8U0BQxNR04KmWQEQk/xQQFDE1HTRRLYGISH4pIChyajrwVEsgIpJfCghKgKZI9lRL\nICKSPwoISkilD2usWgIRkfxRQFACwnkOQLUEmvRIRCQ/FBCUiJqapsmPNEWymg1ERHJNAUEJCSc/\nquSmAzUbiIjkhwKCEhJtOqjksQmUXCgiknsKCEpMTY3GJlAtgYhI7ikgKEEam0C1BCIiuaaAoERV\n+rDGqiUQEcktBQQlKjqscaX2OlAtgYhI7iggKGFhLcHatYUtR6GolkBEJHcUEJQwTX6kWgIRkVxR\nQFDiogmGldh0oFoCEZHcUEBQBsKmg0odsEi1BCIibaeAoAxU+oBF0VqC5asrsJpERCQHFBCUiUof\nsCisJVDTgYhIdhQQlJFKzicIawlERCQ7CgjKTCXnEyiXQEQkewoIykwl5xOox4GISPayCgjM7Ntm\n1mBmW83sSTM7JM3XTTaz7WZ2ezb7lfRE8wkqrekgrCVY8XqFRUMiIm2UcUBgZpOAy4HZwEHAM8C9\nZta7ldftCfwMqLBbVGFU6iiG1b2qWTxxMZNvn8yS5yowu1JEJEvZ1BDMAhY45250zr0AnAlsAaal\neoGZVQE3ARcCqsttB5U8iuGkAyY1BgVqOhARSU9GAYGZdQRGAA+Gy5xzDngA+EILL50NrHPOXZ9N\nISU7ldzrYOQevopECYYiIunJtIagN9ABeDu2/G2gX7IXmNkXgVOA0zIunbRZtNdBJQUFSjAUEcnM\nDjnajgGu2UKz7sDvgNOdcxsy3eisWbPo2bNnwrIpU6YwZcqUbMtZccJeB6NH+0d9vV9WCaIJhtW9\nKuRNi0hFWbRoEYsWLUpYtmnTpqy2Zb7GP82VfZPBFuBrzrk7IstvAHo6574aW38YsBLYhg8aoKlW\nYhsw1DnX7OebmQ0H6urq6hg+fHj670ZSWr7cBwSLF/umhEqx5LklTL59MosnLmbSARX0xkWkYq1c\nuZIRI0YAjHDOrUz3dRnVEDjnPjGzOmA8cAeAmVnw9/wkL1kFHBBbdgnQHZgBqIG3nUS7Iq5bB9On\nF7pE7SMMAibfPpmRe4xUTYGISArZNBnMAxYGgcEKfK+DrsANAGZ2I/C6c+4C59zHwD+jLzazjfhc\nxFVtKbhkbtIkHwzMmAF9+lROTUGYYKimAxGR1DIOCJxztwRjDlwE9AWeBo5yzr0TrLIH8Gnuiii5\nNH26DwYmT/YJh5WQTxAdm6B/j/7UDKopdJFERIpOVkmFzrmrgKtSPDeuldeeks0+JXf69/f/3nln\n5TQdhLUEIiKSnOYyqEA1NTB/vm86qJRBi6LdEJevrqD+lyIiaVJAUKGmT6+8QYvCbogam0BEpDkF\nBBWs0gYtCmsJQJMfiYjEKSCoYNGpkkePhoYK+NFcM6hGkx+JiCSRq5EKpUTV1DSNZLhiRWX0OoiO\nTaBeByIinmoIJGHQoiuvLHRp2kfY60D5BCIingICAfwgRWHPg0oICqL5BJoRUUREAYFETJ9eWd0R\nawbVaEZEEZGAAgJJMGGC/3fy5MpIMgzd+eKdhS6CiEhBKSCQBNGeBysqoGdezaAa5h81nxn3zlCv\nAxGpaOplIM1U2syI0w+dTp9ufdTrQEQqmmoIJKlKSzJUrwMRqXQKCCSlaJJhuQcFGsVQRCqdmgyk\nRWFzwYwZiX+Xo+gohus2r2P6oWX8ZkVEYhQQSKuiQUGfPr45oVxNOmAS6zavY8a9M+jTrU/jqIYi\nIuVOTQaSlmh3xHKfCGn6odMbawo0VbKIVAoFBJKW+ERI5R4UKMlQRCqNAgJJWzgREpR/UFDdq5rF\nExcDSjIUkcqggEAyUklBwaQDJqnpQEQqhgICyVglBQXRpgMFBSJSzhQQSFYqJSiIjk+goEBEypkC\nAslaPCgo18mQwlkRQUGBiJQvBQTSJuG8BwB3lvGEgQoKRKTcKSCQNquUeQ/iQYG6I4pIOVFAIDlR\nKfMeRIMCdUcUkXKioYslZypl3oPonAeAhjcWkbKggEByKh4UTJjgRzksN2F3xMm3T6Z/j/7UDKop\ncIlERNpGTQaSc9Hmg8GDy7NLorojiki5UUAgeTF9elPvg3LtkhhPMlzy3JICl0hEJHsKCCRvJk0q\n/y6J0aBg8u2TufLJMs6oFJGyphwCyatJk2DduvLOKagZVEP9jHrufPFOZtzr3+j0Q8s0o1JEypYC\nAsm7aKLhjBl+dMOaMsvBq+5V3RgEKCgQkVKkgEDaxfTp0KcPTJ7scwrKMSgAEoKCPt36qEuiiJQM\n5RBIu5k0KXHugyVlmoM3YegEAE2bLCIlRQGBtKvohEiTJ5fnqIbxLonqfSAipUABgbS7mhqor28a\nq6AcawrU+0BESo0CAimI6mrf4wB8TUG5BgX1M+qZf9R8Ztw7gyufvFITIolI0VJAIAVTXV0ZzQfT\nD53eGBQMnj9YQYGIFCUFBFJQ8eaDK68sz1ENw6AA4M4Xy3SUJhEpaep2KAVXXV0ZYxXExymYMHQC\n1b3KbJQmESlZqiGQohGf/6AcJ0WKNx+oB4KIFAsFBFJUovMflOtYBdMPnc7iif5NqgeCiBQLNRlI\n0ZkUDO43ebJ/rFvX1KRQLiYdMImRe4xMmP9ATQgiUkgKCKQoTZoEI0f6WRLLdWKk+PwHM+6dwbKT\nl1EzqMySJ0SkJKjJQIpWmGwY9kAYPLh8mxDCHgijF47WcMciUhAKCKToRZMNy3W8gmhegYY7FpFC\nUEAgJWHSpObjFZSbSQdMajbcsQYxEpH2ooBASkZ0uONyHcQoPtyxuiaKSHtRQCAlpbo6saagHPMK\nwmTDaNfEOQ/NKWyhRKTsqZeBlJww2bBPn8SuieXWCyHsmrjw7wupfbQWgJMPOlldE0UkL1RDICUr\nnldQrrUFc8bNYfao2dQ+WqsmBBHJm6wCAjP7tpk1mNlWM3vSzA5pYd3TzGy5mb0XPO5vaX2RTIS1\nBeXeC2HOuDnNRjdU90QRyaWMAwIzmwRcDswGDgKeAe41s94pXjIa+D0wBjgUWAPcZ2b9symwSDKT\nJjVNpVyuCYeTDpiUkHA4euFo5RaISM5kU0MwC1jgnLvROfcCcCawBZiWbGXn3H87565xzj3rnHsJ\nOC3Y7/hsCy2STHwq5XJtQph+6HSWnbyssRlhzkNz1D1RRNoso6RCM+sIjAB+HC5zzjkzewD4Qpqb\n6QZ0BN7LZN8i6aiUhMOaQTWNQxzXPlpL7aO1LJ64mJF7jFTSoYhkJdMagt5AB+Dt2PK3gX5pbuMy\n4A3ggQz3LZK2Skg4hOa5BUo6FJFs5arboQGu1ZXMzgNOAEY75z5ubf1Zs2bRs2fPhGVTpkxhypQp\n2ZZTKkiq2oJynTlxxesrmHz7ZCbfPpl1m9cxrN8wTZQkUuYWLVrEokWLEpZt2rQpq22Zc63ex5tW\n9k0GW4CvOefuiCy/AejpnPtqC689B7gAGO+c+3sr+xkO1NXV1TF8+PC0yyeSyvLlMHq0///8+TBs\nmM85KDcNGxoSplSef9R8TassUmFWrlzJiBEjAEY451am+7qMmgycc58AdUQSAs3Mgr//mup1ZnYu\n8APgqNaCAZF8iCccjh5dnj0R4qMcRoc/VuKhiLQkm14G84AzzOwkM/sMcA3QFbgBwMxuNLPGpEMz\n+x4wF98L4TUz6xs8urW59CIZCJsQli1LzC0ox3ELwi6KyfILFBiISDIZ5xA4524Jxhy4COgLPI3/\n5f9OsMoewKeRl5yF71WwNLap2mAbIu2qpsY/wtyCGTPg3Xfh5JPLqydCda9qqntVN8svAN+UMP3Q\nMkumEJE2ySiHoL0oh0DaS0MDLFwItX6qABYvhpEjyyswCDVsaGgMDABmj5rNuMHjlHgoUmbaJYdA\npNxUV8OcOYlDH5drF8XqXtWNTQnhoEajF47WMMgiAmi2QxHAj1swciSsWJHYRbEceyOEEyaNGzyO\nZ956JqFHgroqilQuBQQigepq/xg5Eu680+cWgE9ALLeRDqFptMNh/YYlBAZqShCpTMohEElhyRJf\nUxBavNjXJJSr5auX81D9Q9Q+6hMqwhqDAT0HaBwDkRKSbQ6BaghEUqikZgRoqjHYd7d9Wbd5XWON\nAaB5EkQqgAICkRa01IxQroHBpAN8NciwfsNY+8HahO6KCgxEypcCApE0hIMaDRsGzzxT/vkFQGMO\nwcg9RjYOhxwNDPr36K88A5EyohwCkSwkyy8o1/ELQstXL2+sMQiFeQaAggORIqEcApF2lCy/AHxg\n0L9/eTYlRGsM1mxak9AzAdRtUaTUqYZApI0aGpoCg9Ds2eU3FHIyLdUaqHeCSGGohkCkQKKJh2vW\nwEMP+aGQa2vLvymhtVqDcHKlMFFRRIqXAgKRHAkDg5oa2HVXn3hYCU0J0DSRUjjQUbx3wqp3VjFu\n8DhAuQYixUpNBiJ5UslNCeAnU1qzaU3CYEegRESRfFOTgUiRaa0poZxrDCCx1iCsHYg3KYTDJIOC\nA5FCUw2BSDu68sqmMQygaYCjAQMqo9YAaJxZUTUHIvmhGgKREjB9uh/IaM2axAGOoLwHOYoKb/bh\nMMn9e/RXzYFIEVANgUgBLV8Oa9c2H+So3JsTkkmn5kBdGUVapxoCkRIU3vSjcyWEwUGlNSekU3MA\naloQyRfVEIgUkYaG1M0J5TqZUmvCmoP4AEig2gORZFRDIFIGomMZDBvW1JwQn2WxUmoNILEWYOQe\nI1nx+grVHojkgQICkSIVbU5YsQLWrUusNVjsBwFkUgUNAhh2ZQQaB0GCptqDZHMrhFSDINIyBQQi\nRS6sNYDEWoMw12DVKhg3rrJqDUKZ1B5A07TN8deKiHIIREpSmGsQDnYUCnsoVGJwEBfmHkDzAZGg\nqQZh7Qdr6d+jvwIEKRvKIRCpINFcg3HjmtcaQGXmG0RFb/DR5gVoOUBItQ2RcqeAQKTERXMN1qxp\nnogIqjmA1AHC2g/Wsm7zuqRNDNEBkkB5CFLeFBCIlIlorkGYiAjNaw4qMRkxmfiv/3jtQDhAUnSQ\nJGjKQ1BTg5QbBQQiZSgeHERrDuLJiFDZNQeh+I09OkAS0GxK51CypgbVJEgpUkAgUuaSBQfRmRdD\nYbMCVOYASMlMOiCxGmXkHiNZs2kN0HJTAyT2aAipNkGKmXoZiFSo5UESfnwuBWhKSATVHrQm2psB\nkicshuI5CaDaBMm9bHsZKCAQkcZujMmCA2jKOxg5UsFBOuJBAjSftClq8UR/gKM1CgoUJFvqdigi\nWUvWrADJuzNGmxZUe5BcsqaBeE4CpM5LCMUTGFvavkhbqYZARFqUasKlUFh7UIlTNudCw4YG1mxa\nk3DTTzaRU1Q0kTH6OtUqCKiGQETyJDoI0oQJLdceRHMPQDUI6YjOzxAVT2AMb/ot5ShA6loFUM2C\ntEwBgYikLdq0AIldGuOTL4WUoJidVIFCfNTF8MbfWvMDNB+uOUq1C6KAQESyFg8QorUDyUZMhMQc\nBFAzQzZS/dJPVasArdcsQPLkxnT2K+VBOQQikjcNDX7ExP79U/dgUDND+wl7P8SDhdZyFkLJBmGK\nbk+1DMVB3Q5FpOiFCYqQOkCA5kGCahHyL1lyYyjdgAGaD8ik5on2p4BAREpONECA9GoR1q5Vj4ZC\nCAOGZNZ+sBYg7aABWm6eUI1D26iXgYiUnHgOAjRNzBTmGaTq7hivRQA1N+RTqiTHqGgOQyhV80Qu\naxxCynFoGwUEIlJU4kFCTU3zZMVUPRqgedJiuA3Jv3SCBmgKHFLd3KM1DpnUOiQbGjq6TXXDbJkC\nAhEpeslu6PHagbAmIVWTQ58+ChSKRaaBQ1SqICLVdNUtmX/UfPp065OyxqHS8h8UEIhISYrfzOM1\nCaFUTQ6QOlBQ00NxSDdwgORDQ0dl0w0zlda6Z5Zqs4YCAhEpG8l+8YeBQpiMGGopUIDEIZmjFCwU\nr/h01S0JB3hq6ebd1vyHZFrrutlaOfIZUKiXgYhUrOXLmwcKLXWHDKUKFkDNEOWupe6ZobZ23WxJ\nPKBItq+Xnn+J0//rdFAvAxGR9KS6eUeHZE4WLLQUMKRqhghfPyn9H7FShDJpxohLlhMRSifASLuZ\n482siqeAQEQkLll3yFCqYAFab4YAWLUKxiVJhI9uT7UM5aktwQQ0n8cCWqgh+PXpGW9fAYGISAZa\nChZS5SuEHnoIamv9oyXxMRZSbU+BQ+VJJ4eg+3vds9q2AgIRkRxq6SZdU5O8dgCabvrp1DKEWmue\n0CRSkgkFBCIi7ai1G3Oy7pPZNk9EJRv+ORkNDV25FBCIiBSZdG7GrTVPRJdnGjxA+s0WyZ5TMFGa\nFNPUzg0AAAYxSURBVBBIo0WLFjFlypRCF6Oi6Ji3v3I65uneeJMN/9zSzb2loaHTEa+N+Mc/FvGl\nLzU/5i2VQ0FF+8sqIDCzbwPnAP2AZ4Dpzrm/tbD+14GLgEHAS8B5zrl7stm35E85XShLhY55+6vU\nY57pDTadZotkzyWvjVjE3LmZHfNkk1dlUo5kFGS0LOOAwMwmAZcDZwArgFnAvWa2j3NufZL1vwD8\nHvg+cBdwIvBHMzvIOffPthReRETyI9ubZ7LaiPnz4Sc/ab5usht4ONZDW2ooUkkWZGQbYLT2XCmO\nN5FNDcEsYIFz7kYAMzsT+DIwDfhpkvXPBu5xzs0L/p5tZkcC3wG+lcX+RUSkiMWDiZtvzizACMd6\nSCabm3Q+g4xUko03ka/gI/7cSy9lVtZQRgGBmXUERgA/Dpc555yZPQB8IcXLvoCvUYi6F/hKJvsW\nEZHK0NJYD9lKFWTk4yad7ngTxSbTGoLeQAfg7djyt4GhKV7TL8X6/VrYT2eAPyz/AyvfbD4M8/rN\n6+ndrXfSF2bzXLFvr732tWnTJlauTHvYa8kBHfP2p2Pe/orlmHdPMl7P3nunXj/b5449FvbYo/ny\n9euhd/JLck6fe/XVVVx8MRDcS9OVq14GBmQyS1Jr6w8CuHjWxW0okmRjxIgRhS5CxdExb3865u1P\nx7wgBgF/TXflTAOC9cA2oG9seR+a1wKE3spwffBNCv8PWA18mGEZRUREKllnfDBwbyYvynj6YzN7\nEnjKOXd28LcBrwHznXM/S7L+YqCLc+4rkWWPA88455RUKCIiUgSyaTKYByw0szqauh12BW4AMLMb\ngdedcxcE6/8CWGZm/4PvdjgFn5iY+VRMIiIikhcZBwTOuVvMrDd+oKG+wNPAUc65d4JV9gA+jaz/\nhJlNAS4JHi8DX9EYBCIiIsUj4yYDERERKT9VhS6AiIiIFJ4CAhERESm+gMDMvm1mDWa21cyeNLND\nCl2mcmZmo8zsDjN7w8y2m9mxhS5TOTOz881shZm9b2Zvm9kfzGyfQpernJnZmWb2jJltCh5/NbOj\nC12uShKc99vNbF7ra0s2zGx2cIyjj4xy9YoqIIhMnDQbOAg/k+K9QRKj5Ec3fGLot8lscCnJzijg\nSuDzwJeAjsB9ZtaloKUqb2vwk6uNCB4PAX8ys30LWqoKEfyoOx1/PZf8eh6f7N8veByeyYuLKqkw\nxRgHa/BjHCSbOElyyMy2A8c55+4odFkqRRDsrgNqnHOPFbo8lcLM3gXOcc5dX+iylDMz6w7UAWcB\nPwL+7pz7n8KWqjyZ2Wx8D77h2W6jaGoIIhMnPRgucz5aaWniJJFStzO+Zua9QhekEphZlZlNxo+d\n8kShy1MBfgX82Tn3UKELUiH2Dpp/XzGzm8xsQCYvztVcBrmQzcRJIiUrqAH7OfCYxuXILzP7LD4A\n6Ax8AHzVOfdCYUtV3oLA60Dg4EKXpUI8CUwFXgT6A3OA5Wb2Wefc5nQ2UEwBQSqZTpwkUiquAvYD\nvljoglSAF4Bh+BqZrwE3mlmNgoL8MLM98MHuEc65TwpdnkrgnIvOW/C8ma0AXgVOANJqGiumgCCb\niZNESpKZ/RI4BhjlnFtb6PKUO+fcp0B98OdKMxsJnI1v25bcGwHsBtQFNWHga4BrzOw7QCdXTAls\nZcg5t8nMXgL2Svc1RZNDEESRdcD4cFlwIo0ng+kbRYpdEAx8BRjrnHut0OWpUFVAp0IXoow9AByA\nbzIYFjz+D7gJGKZgIP+ChM4hQNo/OIqphgBamThJcs/MuuEjyDCKH2xmw4D3nHNrCley8mRmV+En\n+DoW2GxmYY3YJuecpvrOAzO7BLgH32OpB35q9dHAkYUsVzkL2qwT8mLMbDPwrnNuVWFKVd7M7GfA\nn/HNBLsDtfh5hRalu42iCgjSmDhJcu9g4GF8nobDjwMBsBCYVqhClbEz8cf5kdjyU4Ab2700laEv\n/tj2BzYBzwJHKvO93alWIL/2AH4P7Aq8AzwGHOqcezfdDRTVOAQiIiJSGEWTQyAiIiKFo4BARERE\nFBCIiIiIAgIRERFBAYGIiIiggEBERERQQCAiIiIoIBAREREUEIiIiAgKCERERAQFBCIiIgL8f0my\natSkRoqEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb48bd73128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = plt.plot(x, exponential_dist(x,l),'.', markersize = 1)\n",
    "_ = plt.plot(x, pareto_dist(x, xm, alpha), '.', markersize = 1)\n",
    "_ = plt.legend([\"Exponential\", \"Pareto\"])\n",
    "_ = plt.title(\"Comparing the pareto and exponential distributions\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_tail = linspace(0,100,10000)\n",
    "_ = plt.plot(x_tail, exponential_dist(x_tail,l),'.', markersize = 1)\n",
    "_ = plt.plot(x_tail, pareto_dist(x_tail, xm, .01), '.', markersize = 1)\n",
    "_ = plt.legend([\"Exponential, lambda = 1\", \"Pareto, alpha = .01\"])\n",
    "_ = plt.title(\"The heavy tail eventually wins (for all alpha)\")\n",
    "_ = plt.axis([0,100,0,.03])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LogNormal Distribution\n",
    "\n",
    "### A quick note on: the normal distribution\n",
    "\n",
    "It looks like this 😁:  \n",
    "$ f(x) =\\frac{1}{\\sqrt{2\\sigma^2\\pi}} e^{-(x - \\mu)^2 / \\, 2\\sigma^2} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def normal_dist(x, sigma = 1, mu = 0):\n",
    "    '''\n",
    "    x: iterable of x values\n",
    "    sigma: variance\n",
    "    mu: mean\n",
    "    '''\n",
    "    return [1/(sqrt(2*(sigma**2)*pi)) * exp(-(point - mu)**2)/(2*(sigma**2)) for point in x]\n",
    "\n",
    "mu = 0\n",
    "sigma = 1\n",
    "x_mirrored = hstack([-x,x])\n",
    "_ = plt.plot(x_mirrored, normal_dist(x_mirrored, sigma, mu), '.', markersize = 1)\n",
    "_ = plt.ylabel(\"Probabilty mass at x\")\n",
    "_ = plt.xlabel(\"x\")\n",
    "_ = plt.title(\"Normal Distribution\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LogNormal Distribution\n",
    "\n",
    "The LogNormal distribution is, maybe unsurprisingly, a distribution whose logarithm is normally distributed. \n",
    "\n",
    "The PDF of the LogNormal is: $ f(x) = \\frac{1}{x\\sigma\\sqrt{2\\pi}} \\, \\text{exp} [ -\\frac{(\\text{ln}x - \\mu)^2}{2\\sigma^2} ], x>0 $\n",
    "\n",
    "The LogNormal distribution is a heavy tailed distribution, a long tailed distribution, and a subexponential distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def lognormal_dist(x, sigma, mu):\n",
    "    '''\n",
    "    x: iterable of x values\n",
    "    sigma: parameter sigma\n",
    "    mu: parameter mu\n",
    "    '''\n",
    "    pdf = []\n",
    "    for point in x:\n",
    "        pdf.append((1/(point*sigma*sqrt(2*pi))) * exp(-(log(point) - mu)**2 / (2*sigma**2)) )\n",
    "    return pdf\n",
    "mu = 0\n",
    "sigma = 1\n",
    "x_nonzero = x[x != 0]\n",
    "_ = plt.plot(x_nonzero, lognormal_dist(x_nonzero, sigma, mu), '.', markersize = 1)\n",
    "_ = plt.ylabel(\"Probabilty mass at x\")\n",
    "_ = plt.xlabel(\"x\")\n",
    "_ = plt.title(\"LogNormal Distribution\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "_ = plt.semilogx(x_nonzero, lognormal_dist(x_nonzero, sigma, mu), '.', markersize = 1)\n",
    "_ = plt.title(\"'The logarithm of x is normally distributed'\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Central limits\n",
    "\n",
    "The central limit theorem that all undergraduates know and love: the sum of independent, identically distributed random variables goes to the normal distribution.\n",
    "\n",
    "In a later tutorial, I hope to get into when that breaks down theoretically (when those iid random variables have infinite varience, the limit is, in fact, a heavy-tailed distribution and not a normal distribution at all) and when it breaks down practically (obviously, real data cannot have infinite variance, but it can have very, very large variance).\n",
    "\n",
    "For now, I'm going to use the central limit theorem that we know and love to look at the limit of multaplicative (not additive) processes. **Multaplicative processes almost always lead to heavy tails**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# samples of a process\n",
    "sums_of_uniform_samples = []\n",
    "for i in range(0,1000):\n",
    "    sums_of_uniform_samples.append(sum([uniform(1,2) for x in range(0,1000)]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "_ = plt.hist(sums_of_uniform_samples,25)\n",
    "_ = plt.title(\"The Central Limit Theorem works! Hooray!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# now, instead lets multiply together samples of the same (uniform) process\n",
    "multiply_uniform_samples = []\n",
    "for i in range(0,1000):\n",
    "    multiply_uniform_samples.append(reduce(lambda x, y: x*y, [uniform(1,2) for x in range(0,1000)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multaplicative central limit\n",
    "\n",
    "Recall that $ x \\cdot y = \\log(x) + \\log(y) $\n",
    "\n",
    "If a multaplicative process is a process where we are multiplying together many iid random variables, then:\n",
    "\n",
    "$ P_n = X_1 \\cdot X_2 \\cdot X_3 \\cdots \\cdot X_n $  \n",
    "$ \\log(P_n) = \\log(X_1) + \\log(X_2) + \\log(X_3) + \\cdots \\log(X_n) $  \n",
    "And applying the central limit theorm, $\\log(P_n) \\sim N(\\mu, \\sigma^2) $ for some $\\mu, \\sigma^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# note that I'm taking the log of the data before histogramming it \n",
    "_ = plt.hist([log(x) for x in multiply_uniform_samples],50)\n",
    "_ = plt.title(\"The LogNormal Distribution. Whoa.\")\n",
    "_ = plt.xlabel(\"log(x)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting data\n",
    "\n",
    "My last note (and here there will be much more explanation in future notebooks), is on fitting a little bit of real data. Because real data is fun."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# This is a random list of follower counts for people on Twitter\n",
    "followers = []\n",
    "for follower_count in sorted(open(\"random_selection_of_followers.txt\",\"r\").read().split(\"\\n\")):\n",
    "    try:\n",
    "        followers.append(int(follower_count))\n",
    "    except ValueError:\n",
    "        pass\n",
    "# I'm grabbing exacly 1 million of those counts, because 1 million is a nice tidy number\n",
    "followers = sorted(sample(followers, 10**6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common way of plotting these kinds of distributions is a rank-oreder plot, or (pictured below) a complementary cumulated distribution function (CCDF) with the idea that we want to visualize the probability that $X$ has frequency greater or equal to $x$.\n",
    "\n",
    "I simply order the observations of the number of followers and plot the fraction of users greater than $x$ against $x$ (the actual follower count)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "_ = plt.loglog(followers, [1 - i/10**6 for i in range(0,10**6)], '.', markersize = 1)\n",
    "_ = plt.title(\"Follower counts on Twitter\")\n",
    "_ = plt.ylabel(\"Probabilty that the number of followers \\n of a given user is greater than x\")\n",
    "_ = plt.xlabel(\"x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# since follow counts are discrete, it's easy to bin them by making each bin 1 follower long\n",
    "freq_of_follow_counts = sorted(Counter(followers).items(), key = lambda x:x[0])\n",
    "# follow counts\n",
    "follow_count = [x[0] for x in freq_of_follow_counts]\n",
    "# frequency of each follower count\n",
    "frequency = [x[1] for x in freq_of_follow_counts]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "_ = plt.loglog(follow_count, frequency, '.')\n",
    "_ = plt.title(\"# of followers of random Twitter users\")\n",
    "_ = plt.xlabel(\"# of followers x\")\n",
    "_ = plt.ylabel(\"# of people with x followers\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The powerlaw Python library\n",
    "\n",
    "A nice person made this thing. It's not the final word on fitting long tailed distributions to empirical data, but it's based on solid statistics and it enables us to fit power law and lognormal distributions to the same dataset and compare them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# run a \"fit\" on the vector of all 1 million follower numbers\n",
    "results = powerlaw.Fit(followers)\n",
    "# https://github.com/jeffalstott/powerlaw/issues/28 <- I know you get a warning, \n",
    "#                                                       he knows you get a warning, it's not a big deal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"The alpha value of the Pareto distribution that fits this data best: {}.\".format(results.power_law.alpha))\n",
    "print(\"The x_m value of the Pareto distribution that fits this data best: {}.\".format(results.power_law.xmin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"The sigma value of the LogNormal distribution that fits this data best: {}.\".format(results.lognormal.sigma))\n",
    "print(\"The mu value of the LogNormal distribution that fits this data best: {}.\".format(results.lognormal.mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Importantly, this package allows you to compare the powerlaw and lognormal distributions\n",
    "'''\n",
    "R is the log-likelihood ratio between both distributions and p is the p-value. \n",
    "Generally, we can say that if R is positive, than the first distribution is a better fit while if it is negative, \n",
    "the second one should be prefered. \n",
    "'''\n",
    "R, p = results.distribution_compare('power_law', 'lognormal') \n",
    "\n",
    "print(\"The R value: {}\".format(R))\n",
    "print(\"The p-value: {}\".format(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Those results don't look very good, but sometimes real distributions can seem mroe piecewise\n",
    "# let's try fitting jsut the top piece of that distribution\n",
    "results_without_tail = powerlaw.Fit([i for i in followers if i <= 10**2.5]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "R_without_tail, p_without_tail = results_without_tail.distribution_compare('power_law', 'lognormal') \n",
    "\n",
    "print(\"The R value: {}\".format(R_without_tail))\n",
    "print(\"The p-value: {}\".format(p_without_tail))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"The sigma value of the LogNormal distribution that fits this data best: {}.\".format(results_without_tail.lognormal.sigma))\n",
    "print(\"The mu value of the LogNormal distribution that fits this data best: {}.\".format(results_without_tail.lognormal.mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cuttoff_index_number = [ n for n,i in enumerate(followers) if i > 10**2.5 ][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "results_1.lognormal.plot_pdf()\n",
    "plt.loglog([i for i in follow_count], [i/10**6 for i in frequency], 'r.', markersize = 1)\n",
    "plt.title(\"LogNormal fit on empirical data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "results_1.power_law.plot_pdf()\n",
    "plt.loglog([i for i in follow_count], [i/10**6 for i in frequency], 'r.', markersize = 1)\n",
    "plt.title(\"Powerlaw fit on empirical data\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
